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I.  INTRODUCTION  AND  SUMMARY 

1.1  INTRODUCTION 

The  objectives  of  this  contract  may  be  summarized  as 
follows : 

1 . Determine  a source  model  for  the  1975  Pocatello 
Valley  earthquake. 

2.  Compute  the  likely  ground  motion  at  the  Wing  V 
Minuteman  site  due  to  the  Pocatello  Valley  event 
and  determine  those  features  of  the  source  or 
propagation  path  that  may  have  caused  this  ground 
motion  to  be  peculiar  or  untypical  of  western 

U.  S.  earthquakes. 

3.  Using  synthetic  seismograms  for  typical  western 
U.  S.  earthquakes  and  earth  structure,  determine 
the  important  characteristics  of  the  ground  motions 
at  regional  distances  within  200  km  of  these  events. 

4.  Provide  initial  estimates  of  the  similarities  and 
differences  between  the  strong  ground  motions 
generated  by  earthquakes  and  large  surface  nuclear 
■z^xplosions. 

In  this  report  we  will  be  concerned  with  the  first 
objective  and  that  portion  of  the  second  objective  which 
compares  the  source  model  of  the  Pocatello  Valley  earthquake 
with  other  western  U.  S.  earthquakes.  A companion  report 
(Rodi,et  a^. , 1979)  addresses  the  remaining  objectives  with 
special  emphasis  on  the  ground  motion  at  Wing  V from  both 
the  Pocatello  Valley  earthquake  and  the  1975  Yellowstone 
earthquake . 

A variety  of  techniques  were  used  in  order  to  carefully 
delineate  the  source  mechanism  for  the  initial  faulting  of  the 
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1975  Pocatello  Valley  earthquake.  The  data  include  after- 
shock locations  and  observations  of  long  period  and  short 
period  P waves.  We  use  analytical  models  and  a three-dimen- 
sional finite  difference  earthquake  model  to  compute  synthetic 
seismogrcuns  to  match  these  data.  We  find  that  the  Pocatello 
Valley  earthquake  is  characterized  by  normal  faulting  on  a 
plane  dipping  39°  to  the  west  with  a large  left-lateral 
strike  slip  component.  This  event  turns  out  to  be  quite 
similar  in  most  respects  to  the  1971  San  Fernando  earthquake 
in  Southern  California. 

1.2  SUMMARY 

The  technical  discussion  in  this  report  is  divided  in- 
to three  sections.  In  Section  II  we  use  analytical  models  to 
infer  a fault  model  for  the  Pocatello  Valley  earthquake.  As 
part  of  this  contract/  ILLIAC  IV  computer  time  was  provided 
in  order  that  a three-dimensional  finite  difference  model  of 
faulting  could  be  exercised.  Section  III  is  devoted  to  a dis- 
cussion of  that  work  and  its  implications  for  the  physics  of 
earthquake  faulting.  In  Section  IV  the  results  of  the  analy- 
tical modeling  of  the  Pocatello  Valley  earthquake  and  the 
finite  difference  calculations  are  compared  and  combined  to 
complete  our  source  model  for  the  Pocatello  Valley  event. 
Important  conclusions  about  the  nature  of  earthquake  fault- 
ing in  general  are  also  drawn. 

The  techniques  applied  in  Section  II  to  infer  the 
source  parameters  of  the  Pocatello  Valley  earthquake  include 
the  following: 

• A focal  mechanisms  solution  is  constructed  from 
observations  of  P wave  first  motion. 

• Results  from  a study  of  the  aftershock  sequence 
by  Arabasz,  et  al.  (1975b)  were  examined  for 
indications  of  the  fault  plane  orientation. 
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• Using  the  Archambeau/Minster  earthquake  source 
model  (Archambeau,  1968;  Minster,  1973)  and 
methods  for  computing  synthetic  seismograms  in 


layered  earth  models  (Bache  and  Harkrider,  1976) , 
synthetic  seismograms  were  computed  for  compari- 
son to  long  and  short  period  P wave  observations. 

• The  model  parameters  were  varied  to  obtain  esti- 
mates for  the  source  dimension,  rupture  velocity 
and  stress  drop  as  well  as  confirming  the  orienta- 


tion. The  main  diagnostic  features  of  the  obser- 
vations are  the  waveform  and  frequency  content  of 
the  P and  pP  phases,  their  time  delay  and  their 
relative  amplitudes. 

An  important  conclusion  of  Section  III  is  that  no  uni- 
form stress  drop,  uniform  rupture  velocity  model  can  simul- 
taneously fit  the  long  and  short  period  data.  If  the  long 
period  level  (moment)  is  large  enough,  the  synthetic  short 
period  records  are  much  too  large  — by  a factor  of  3 to  5. 

We  were  led  to  the  same  conclusion  in  an  earlier  study  of  the 
1971  San  Fernando  earthquake  (Bache  and  Barker,  1978) . 

To  improve  the  fit  to  the  data  over  a broad  frequency 
range,  variable  stress  drop  and  variable  rupture  velocity 
effects  were  introduced  into  the  analytical  source  model. 

This  is  done  in  a rather  ad  hoc  fashion  and  the  results  are 
far  from  unique.  However,  we  claim  that  they  are  qualita- 
tively correct.  The  new  source  model  gives  a reasonably  good 
fit  to  the  long  and  short  period  data  for  the  first  5 to  8 
seconds  of  faulting. 

The  approach  in  Section  III  is  entirely  different.  In 
this  section  our  main  attention  is  to  the  physics  of  earth- 
quake faulting  for  events  like  the  Pocatello  Valley  earth- 
quake. Two  three-dimensional  finite  difference  calculations 
were  done  to  simulate  faulting  in  a uniform  whole  space.  The 

I calculations  were  done  on  the  ILLIAC  IV  computer  using  the 

/ * 

I TRES  computer  program  (Cherry,  1977).  The  first  calculation 

was  for  a linearly  elastic  medium  while  the  second  allowed 
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elastic-plastic  material  behavior  near  the  fault  plane. 

Using  the  method  of  Bache  and  Harkrider  (1976)  , the  results 
of  both  calculations  were  expressed  in  a form  compatible  with 
the  progreutis  for  computing  synthetic  seismograms.  The  far- 
field  displacements  and  displacement  spectra  were  studied. 

The  accuracy  of  the  finite  difference  method  was  verified 
by  comparing  the  results  to  analytical  and  numerical  solutions 
of  crack  problems.  We  find  that  main  effect  of  the  simple 
plasticity  admitted  to  the  second  calculation  was  less 
abrupt  stopping  of  the  rupture  at  the  edge  of  the  fault  plane. 
This  reduces  the  radiation  of  high  frequency  energy.  Plastic 
flow  outside  the  fault  plane  increases  the  seismic  moment. 

In  Section  II  we  develop  a model  for  the  Pocatello 
Valley  earthquake  in  terms  of  an  analytical  model  for  which 
parameter  variations  are  convenient.  However,  interpreta- 
tion of  the  derived  parameters  in  terms  of  earthquake  physics 
is  difficult  for  any  analytical  model  and  is  especially  so 
for  the  Archambeau/Minster  model  which  is  cast  in  terms  of 
an  unrealistic  spherical  geometry.  On  the  other  hcind,  the 
ILLIAC  calculations  of  Section  III  are  close  to  the  most 
detailed  and  realistic  models  of  earthquake  faulting  that 
are  currently  available.  Taken  together,  the  two  can  be 
used  to  define  the  source  with  considerable  confidence.  This 
is  our  objective  in  Section  IV. 

Comparison  of  the  finite  difference  and  analytical 
models  by  means  of  synthetic  seismograms  appropriate  for 
stations  observing  the  Pocatello  Valley  event  support  our 
conclusions  about  earthquake  faulting  reached  with  the 
analytical  models.  That  is,  we  cannot  simultaneously  match 
the  long  and  short  period  data  with  a single  rupture  veloc- 
ity/single stress  drop  source  model. 

Another  interesting  facet  of  the  comparison  is  that 
source  directivity  effects  can  be  seen  in  that  data.  The 
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finite  difference  simulations  are  bilateral  and  give  poorer 
agreement  with  the  short  period  data  than  the  unilateral 
analytical  model. 

Finally,  we  conclude  that  the  variable  rupture  veloc- 
ity analytical  model  derived  in  Section  II  is  a good  model 
for  the  initial  few  seconds  of  rupture  of  the  Pocatello 
Vally  earthquake.  The  long  and  short  period  teleseismic 
records  show  that  the  source  is  more  complex  and  has  longer 
duration.  The  later  faulting  is,  however,  considerably 
smaller  as  a source  of  seismic  waves  so  our  Model  II  should 
give  a reasonable  estimate  for  the  peak  ground  motions  at 
the  regional  ranges  of  ultimate  interest. 
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II.  AN  ANALYTICAL  SOURCE  MODEL  FOR  THE  POCATELLO 


VALLEY  EARTHQUAKE 

2 . 1 INTRODUCTION 

The  March  28,  1975  Pocatello  Valley  earthquaUce  was  an 
6.1  event  occurring  near  the  Idaho-UtaUi  border.  Our  pro- 
cedure for  deducing  the  ground  motions  of  interest  involves 
two  steps.  First,  we  use  availcQsle  data  to  infer  a detailed 
source  model  for  the  event.  Then  we  use  this  source  model  to- 
gether with  estimates  for  the  earth  structure  to  theoretically 
compute  the  ground  motions  of  interest.  This  section  de- 
scribes our  first  step;  inference  of  the  source  from  analyt- 
ical models. 

Previous  work  on  the  Pocatello  Valley  event  has 
included  a thorough  study  of  the  aftershock  sequence  (Arad^asz, 
et  al. , 1975b) . Fault  plane  solutions  were  constructed  by 
A.  Rogers  (USGS,  private  communication)  from  the  first  motion 
P wave  observations  and  by  Battis  and  Hill  (1977)  from  far- 
field  surface  wave  spectra.  We  take  advantage  of  this  work 
in  our  more  detailed  study. 


The  details  of  our  source  model  are  largely  determined 
by  a comparison  of  synthetic  and  observed  short  period  seis- 
mograir-  at  teleseismic  distances;  that  is,  at  ranges  greater 
than  30’.  These  data  constrain  the  source  in  the  frequency 
range  from  0.5  to  2.0  Hz.  Of  course,  the  accuracy  with 
which  we  can  constrain  the  source  in  this  frequency  band  de- 
pends on  our  ability  to  account  for  path  effects.  Unfor- 
tunately, there  are  only  a few  good  teleseismic  recordings 
of  this  event  and  this  increases  the  range  of  uncertainty. 


The  source  model  we  use  for  our  synthetic  seismogram 
computations  is  the  analytical  relaxation  model  of  Archambeau 
(1968)  and  Minster  (1973).  This  model  has  previously  been 
used  to  model  the  1971  San  Fernando  earthquake  (Bache  aind 
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Barker,  1978) . The  analytical  form  is  convenient  for  varying 
the  source  parameters.  In  Section  III  we  will  describe  a de- 
tailed finite  difference  fault  model  for  earthquake  faulting 
that  reproduces  the  fault  physics  with  much  greater  accuracy. 
Then  in  Section  IV  we  will  compare  the  results  of  the  finite 
difference  calculations  with  our  analytical  model  and  discuss 
the  implications  for  the  Pocatello  Valley  event  and  other 
western  U.S.  earthquakes. 

2.2  FOCAL  MECHANISM  DETERMINATION  FROM  P WAVE  POLARITY 

Within  weeks  of  the  event,  Rogers  (1975)  determined  a 
focal  mechanism  solution  from  immediately  available  P wave 
first-motion  readings.  Another  solution  was  obtained  by 
Battis  and  Hill  (1977)  who  used  surface  wave  spectra  to  esti- 
mate the  fault  plane  orientation.  Arabasz,  et  al . (1975b) 

reported  on  a detailed  study  of  the  aftershock  sequence  and 
used  this  information  to  deduce  the  fault  parameters.  These 
three  solutions  are  summarized  in  Table  1. 

Looking  at  the  previous  focal  mechanism  solutions  and 
how  they  were  obtained,  we  decided  to  verify  them  by  repeating 
the  P wave  first-motion  solution,  adding  data  unavailable  to 
Rogers.  For  this  purpose,  first-motion  readings  were  obtained 
from  the  U.S.  Geological  Survey  (USGS)  (Rogers,  private  com- 
munication) and  from  the  California  Institute  of  Technology 
(Stewart,  private  communication) , who  also  computed  the  lower 
hemisphere  location  for  the  ray  to  each  station.  A summary 
of  station  location,  first  motion,  amplitude  and  magnitude 
data  is  given  in  Table  2. 

The  P wave  polarity  data  are  plotted  in  a lower  hemi- 
sphere display  in  Figure  1.  A rather  well-constrained  focal 
mechanism  is  determined  from  the  data.  The  orientation  of 
the  two  orthogonal  best-fitting  solutions  is  indicated  on  the 
figure.  One  of  these  solutions  is  for  normal  faulting  on  a 


TABLE  2 

POCATELLO  EARTHQUAKE,  MARCH  28,  1975 


074) 


POCATELLO  EARTHQUAKE,  MARCH  28,  1975  (continued) 
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o DILITATIONAL  FIRST  MOTION 

• COMPRESSIONAL  FIRST  MOTION 

Solution  1 Solution  2 

Strike  * N1*E  Strike  = N45''E 
Dip  = 60®E  Dip  * 39“W 

Slip  » Up  64  »N  Slip  = Up  53  ®N 

Normal  Faulting  Normal  Faulting 

Figure  1.  The  lower  hemisphere  fault  plane  solution  is  shown 
for  the  Pocatello  Valley,  Idaho  earthquake  of 
28  March  1975. 


fault  plane  dipping  60®  to  the  east  with  some  right-lateral 
strike-slip  component.  The  second  solution  is  for  normal 
faulting  on  a plane  dipping  39°  to  the  west  with  a large 
left-lateral  strike-slip  component.  Choice  between  the  two 
must  be  based  upon  other  information. 

2 . 3 GEOLOGICAL  EVIDENCE 

Arabasz,  ^ al.  (1975b)  reported  on  a detailed  study 
of  the  aftershocks  associated  with  this  event.  More  than 
400  aftershocks  were  accurately  located  and  are  mapped  in 
Figures  2 through  5.  The  focal  mechanisms  are  predominantly 
for  normal  faulting  with  some  strike-slip  component.  The 
depth  distribution  of  the  aftershocks  is  indicated  in  the 
cross-sections  plotted  in  Figures  3 through  5. 

The  conventional  wisdom  in  seismology  is  that  the 
aftershock  distribution  can  be  used  to  indicate  the  location 
of  the  fault  plane.  We  have  plotted  the  main  shock  hypo- 
center  in  Figures  3 and  4 using  the  epicenter  location  of 
Arabasz,  et  al . (1975b)  from  Figure  2.  The  depth  of  8.7  km 
was  deduced  in  the  synthetic  seismogram  studies  to  be  des- 
cribed in  later  sections.  It  is  not  greatly  different  from 
the  preliminary  depth  determination  of  5.0  1cm  given  by 
Arab,  sz,  et  al . (1975a).  There  may  be  errors  in  both  the 

depth  and  epicenter  determinations,  but  the  aftershock  dis- 
tributions on  Figures  2-4  certainly  support  the  chosen 
locations.  Section  2 in  Figure  2 turns  out  to  be  exactly 
aligned  with  the  westerly  dipping  focal  mechanism  solution 
we  have  called  Solution  2 in  Figure  1.  The  depth  distribution 
of  aftershocks  in  this  section  (Figure  4)  is  remarkably  con- 
sistent with  the  39°w  dip  of  the  focal  mechanism  solution. 

The  aftershock  distribution  can  also  be  used  to  help  define 
the  fault  dimensions. 
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Figure  2.  The  locations  are  plotted  for  the  Pocatello  Valley 
earthquake  and  more  than  400  fore  and  aftershocks 
(reproduced  from  Arabasz,  et  al. , 1975b).  The 
strike  for  Solution  2 from  Figure  1 is  indicated. 
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Figure  4.  The  depth  distribution  of  events  is  shown  for 

Section  2.  The  main  shock  location  is  clotted  as 
in  Figure  3 and  the  dip  solutions  of  Figure  1 are 
indicated  by  a solid  line  (reproduced  from  Arabasz 
et  al. , 1975b) . 
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Except  for  finite  source  directionality  effects,  there 
are,  of  course,  no  differences  in  the  seismic  radiation  from 
the  complementary  solutions.  Battis  and  Hill  (1977)  inverted 
surface  wave  spectra  to  find  complementary  solutions  quite 
close  to  ours.  They  state  a preference  for  a solution  dipping 
60°  to  the  east  because  of  its  consistency  with  past  seismic- 
ity in  the  area  and  particularly  the  1962  Cache  Creek  earth- 
quake. But  this  seems  a very  weaik  argument  for  preferring 
the  east  dipping  solution  over  that  dipping  39°  to  the  west. 

In  summary,  the  aftershock  distributions  strongly  sup- 
port the  fault  plcine  solution  with  a 39 °W  dip.  We  will  later 
show  results  from  our  attempts  to  model  the  short  and  long 
period  teleseismic  body  waves  using  this  orientation. 

2.4  CALCULATION  OF  SYNTHETIC  SEISMOGRAMS 

If  we  are  to  use  far-field  ground  motions  to  deduce 
the  important  characteristics  of  the  earthquake  source,  we 
must  have  a source  model  and  an  accurate  means  for  accounting 
for  the  source-receiver  travel  path.  We  briefly  describe  the 
the  numerical  methods  and  models  used.  For  a more  detailed 
account,  see  Bache  and  Barker  (1978) . 

The  source  model  we  are  using  for  this  event  is  the 
Archambeau/Minster  model  (Archambeau,  1968;  Minster,  1973). 

The  far-field  P wave  spectrum  for  this  model  has  roughly  the 
form  indicated  in  Figure  6.  The  long  period  level  (u^)  and 
corner  frequency  (f  ) are  proportional  to  stress  drop  (Ao) , 
fault  dimension  (L) , rupture  velocity  (V^)  and  P wave  veloc- 
ity (a^) , as  indicated.  They  are  also  radiation  pattern 
dependent,  though  much  more  weakly.  The  source  model  also 
includes  directionality  effects;  that  is,  enhanced  high-fre- 
quency energy  in  the  direction  of  rupture  propagation. 
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Figure  6 . 
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For  determining  the  amplitude  of  the  source  spectrum 
we  must  account  for  elastic  and  nonelastic  attenuation  in  the 
earth.  Synthetic  seismograms  are  computed  with  the  method 
of  Bache  and  Harkrider  (1976)  for  embedding  a complex  source 
representation  in  a plane-layered  earth  model.  The  source 
region  crustal  model  is  tabulated  in  Table  3 (Keller,  et  al. , 
1975) , while  that  for  the  rest  of  the  path  is  Model  C2  of 
Anderson  and  Hart  (1976) . Another  important  parameter  is  the 
effective  Q to  account  for  anelastic  attenuation  along  the 
path.  The  short  period  waves  are  particularly  sensitive  to 
this  parameter.  We  approximate  the  effect  of  Q by  applying 
the  operator  (Strick,  1970) 

exp  j-  TTf  t*  [l  - I i (1) 

where  f is  frequency  and  t*  is  the  ratio  of  the  travel  time 
to  the  effective  path  Q. 

The  procedure  for  deducing  the  source  parameters  from 
the  far- field  body  wave  recordings  includes  several  steps. 
Given  the  source  orientation,  the  P and  S velocities  in  the 
source  region  and  a t*  for  the  travel  path,  we  guess  values 
for  the  other  source  parameters.  Synthetic  seismograms  are 
then  computed  and  compared  to  the  observations.  The  free 
source  parameters  (Aa , L,  depth)  are  then  adjusted  to 

optimize  the  agreement. 

Because  of  the  nature  of  the  source  model,  the  param- 
eters that  are  actually  fixed  by  the  procedure  are  Poisson's 
ratio  (v),  Aa  and  the  ratios  L/V  , V /B . That  is,  events 

K K 

with  the  same  v,  V^^/S  and  Aa  are  indistinguishable. 

In  the  sequel  we  will  assume  \j  is  known. 

2.5  COMPARISON  OF  SYNTHETIC  AND  OBSERVED  SEISMOGRAMS 

Comparing  synthetic  and  observed  seismograms,  we  have 
constructed  a model  of  the  Pocatello  Valley  earthquake.  The 
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NORTHERN  WASATCH  (Keller,  et  al.,  1975) 


Layer 

Depth 

Thickness 

a 

s 

®s 

"s 

1 

1.4 

1.4 

3.4 

2.0 

2.25 

2 

15.5 

14.1 

5.9 

3.4 

I 

2.73 

3 

20.0 

4.5 

6.0 

3.5 

2.80 

4 

25.0 

5.0 

6.0 

3.5 

2.80 
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event  of  this  size  the  source  corner  frequency  is  well  out- 
side the  bandpass  of  the  long  period  instrument.  That  is, 


these  data  are  mainly  sensitive  to  the  source  orientation, 
depth  and  moment.  Finite  source  effects  such  as  the  source 
dimension  and  rupture  velocity  do  not  play  a role. 

We  begin  with  a simple  model  including  constant  rup- 
ture velocity  and  stress  drop.  This  model  (referred  to  as 
Model  I)  is  characterized  by  the  following  parameters. 

MODEL  I 

Fault  Length,  L:  3 km 

Rupture  Velocity,  V : 3 km/sec 

R 

(Unilateral  fault  propagation  toward  the  surface) 

Stress  Drop,  Aa : 257  bars 

Orientation:  Strike  N45‘’ 

Dip  39 "W 
Slip  Up  53°N 

Of  course,  L and  V have  no  significance  other  than 

R 

the  fact  that  they  give  the  correct  moment.  The  stress  drop 
is  selected  to  bring  the  amplitudes  of  observed  and  synthetic 
records  into  agreement,  as  we  shall  see.  For  the  Archambeau/ 
Minster  model  the  moment  is  (Minster  and  Suteau,  1977) 

M^  = 1.024  L^  Aa  . (2) 

The  observed  and  synthetic  long  period  body  wave 
seismograms  are  compared  in  Figure  7 for  four  WWSSN  stations 
and  the  Canadian  seismic  network  station  RES.  The  agreement 
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is  good  for  the  first  5 to  8 seconds  in  which  the  main  arriv- 
als are  the  direct  P wave  from  the  fault  initiation  and  the 
associated  free  surface  reflections. 

In  Table  4 we  compare  the  observed  and  synthetic 
amplitudes  and  periods  for  the  maximum  peak-to-peak  excur- 
sion. The  period  is  twice  the  time  lag  between  peaks.  This 
comparison  is  done  for  the  five  stations  plotted  in  Figure 
7 and  several  others.  The  period  is  very  difficult  to  measure 
on  the  observed  records  and  the  Tq  in  the  table  may  be  in 
error  by  a second  or  so.  Recognizing  this,  the  observed  and 
synthetic  periods  are  in  quite  close  agreement. 

The  amplitude  ratios  in  Table  4 include  the  effect 
of  differing  instrument  amplification  at  the  periods  Tq  and 
Tg.  If  we  ignore  the  apparent  period  differences  and  assume 
the  instrument  amplification  is  the  same  for  Aq  and  Ag,  the 
mean  is  increased  by  10  percent  and  the  standard  deviation  is 
slightly  reduced.  These  differences  are  not  significant. 

If  the  synthetic  seismograms  are  computed  with  some 
standard  value  for  Mq,  the  inferred  moment  is  then  the  pro- 
duct of  Ao/Ag  and  this  standard  value.  The  inferred  moment 
values  are  also  listed  in  Table  4.  These  values  are,  of 
course,  dependent  on  the  elastic  and  anelastic  properties 
of  the  path  model  chosen.  In  this  case  we  are  using  the 
model  C2  of  Anderson  and  Hart  (1976)  and  t*.  = 0.8.  This 
value  is  near  the  upper  limit  of  values  inferred  from  spec- 
tral analysis  of  Basin  and  Range  events  (e.g.,  Der  and 
McElfresh) , but  near  the  lower  limit  of  values  obtained  by 
ananalysis  employing  synthetic  seismograms  (e.g.,  Bache, 
et  al . , 1975;  Burdick,  1978).  Thus,  it  seems  a reasonable 
choice. 


TABLE 


The  moment  values  inferred  from  the  nine  stations  in 
Table  4 are  remarkably  consistent;  the  standard  deviation  is 
only  28  percent  of  the  mean,  or  less  if  we  ignore  the  fre- 
quency dependent  instrument  correction.  The  mean  value  is 
2 4 

7.1  X 10  dyne-cm.  This  is  near  the  preliminary  moment  esti- 
mate given  by  Arabasz,  et  al.  (1975a)  and  is  almost  the  same 
as  the  value  obtained  by  Battis  and  Hill  (1977)  from  surface 
waves  (see  Table  1) . The  stress  drop  for  Model  I is  Aa  = 253 
bars,  which  scales  the  moment  of  this  value. 

Now  let  us  direct  our  attention  to  the  short  period 
records.  These  data  are  sensitive  to  rather  fine  details 
of  the  source.  Since  the  teleseismic  short  period  records 
are  often  quite  complicated  and  are  influenced  by  both  path 
and  source  details,  they  have  not  often  been  used  to  infer 
source  characteristics.  However,  we  were  successful  in 
modeling  the  short  period  records  for  the  San  Fernando  earth- 
quake CBache  and  Barker,  1978)  and  were  anxious  to  see  how 
well  other  events  could  be  modeled. 

Synthetic  short  period  records  from  Model  I are  com- 
pared to  the  observed  records  at  five  teleseismic  stations 
in  Figure  8.  Synthetic  seismograms  also  are  shown  for  a 
model  we  call  Model  II,  but  we  will  discuss  these  a bit  later. 
There  are  few  stations  where  the  records  are  clear  enough 
to  digitize  and  those  shown  are  the  best  of  them.  The  main 
arrivals  on  the  synthetic  are  P and  pP.  The  sP  phase  is 
rather  small  for  this  source  orientation. 

Comparison  of  the  synthetic  and  observed  records  allows 
us  to  identify  the  largest  phase  as  pP.  The  P wave  is  small 
and  characterized  by  low  signal-to-noise . It  is  not  easy 
to  identify  and  measure  on  the  observations . 

The  main  diagnostic  features  of  the  comparison  in- 
clude the  waveform  and  frequency  content  of  P and  pP  and  the 


Fi<jure  8.  Observed  and  synthetic  seismograms  are  compared  at  five  WWSSN 
stations  for  Models  [ and  II.  Tick  marks  are  at  1 second 
intervals . 


relative  amplitude  of  the  two  phases.  These  contain  infor- 
mation about  the  fault  propagation  in  the  first  few  seconds 
of  rupture  and  the  source  radiation  pattern  (orientation) . 

The  spacing  between  P and  pP  can  be  used  to  fix  the  depth. 
Precision  in  this  depth  estimate  depends  on  our  matching  the 
waveform  well  enough  to  unambiguously  pick  the  onset  of  the 
two  phases. 

Comparing  the  observed  waveforms  to  those  from  the 
Model  I synthetic,  we  see  that  the  two  are  in  reasonable 
agreement,  especially  when  we  take  the  rather  poor  quality 
of  the  data  into  account.  The  shape  of  the  pP  phase  (which 
also  includes  a small  contribution  from  sP)  is  roughly  cor- 
rect. The  observations  are,  perhaps,  a bit  longer  period 
than  the  synthetics.  It  is  difficult  to  judge  how  well  we 
have  fit  the  waveform  of  the  P wave  since  it  is  not  very 

clear  on  the  records.  The  spacing  between  P and  pP  seems  to 

be  fit  about  as  well  as  possible.  The  relative  amplitude  of 

the  P and  pP  phases  are  fit  fairly  well,  though  the  P phase 

is  a bit  too  large  at  several  stations. 

In  Table  5 we  compare  the  amplitudes  of  the  observed 
and  synthetic  records.  The  source  is  Model  I with  Aa  = 253 
bars;  that  is,  the  source  that  gave  the  best  fit  to  the  long 
period  data.  If  our  source  fit  the  long  and  short  period 
data  equally  well,  the  amplitude  ratios  would  have  a mean  of 
uni  ty . 

In  Table  5 we  list  observed/synthetic  amplitude  ratios 
for  both  the  P phase  and  the  maximum  phase  that  we  have 
identified  as  being  primarily  pP.  We  have  attempted  to  com- 
pare amplitude  measurements  made  from  corresponding  portions 
of  the  waveforms  and  where  the  measured  periods  are  approxi- 
mately the  same.  Small  bars  on  the  synthetic  records  indi- 
cate where  the  peak-to-peak  amplitudes  and  periods  were 
measured.  It  is  important  to  compare  amplitudes  where  the 


apparent  frequency  content  is  about  the  same  to  avoid  having 
to  make  a large  correction  for  the  frequency-dependent 
instrument  response.  This  instrument  correction  is  a rough 
approximation  at  best,  and  we  try  to  minimize  its  effect. 

The  P wave  amplitude  ratios  in  Table  4 have  a mean  of 
0.16  while  the  amplitude  ratios  for  the  maximum  amplitude 
have  a mean  of  0.25.  For  the  former  the  scatter  is  very 
large  with  the  standard  deviation  being  83  percent  of  the 
mean  for  this  small  population.  The  maximum  phase  data  have 
much  less  scatter,  which  is  probably  due  to  this  phase  being 
much  better  defined. 

The  data  in  Table  4 force  us  to  the  conclusion  that 
Model  I cannot  simultaneously  fit  the  long  and  short  period 
data.  When  we  adjust  the  model  to  have  the  correct  moment, 
the  synthetic  short  period  body  waves  are  too  large  by  a 
factor  of  4 to  6.  Yet  we  have  matched  the  frequency  content 
of  the  short  period  records  reasonably  well. 

The  only  difference  between  the  long  and  short  period 
synthetic  records  is  in  the  instrument  response.  Both  use 
the  same  source  and  earth  models.  The  discrepancy  between 
short  and  long  period  sources  could  be  reduced  a bit  by 
supposing  a large  t*,  but  this  could  only  account  for  a 
relatively  minor  adjustment. 

To  match  the  data,  we  need  a source  that  has  the  same 
long  period  level  as  Model  I,  but  4 to  6 times  smaller  ampli- 
tudes in  the  0.5  to  2.0  Hz  range  sampled  by  the  short  period 
instrument.  Since  the  frequency  content  of  Model  I is  about 
right,  we  cannot  do  this  by  merely  supposing  a larger  source 
dimension  (with  lower  stress  drop)  to  reduce  the  corner 
frequency.  ^ 

In  our  study  of  the  San  Fernando  earthquake  (Bache  and 
Barker,  1978)  we  found  the  same  kind  of  difficulty  in 
simultaneously  matching  long  and  short  period  data.  We 


concluded  that  this  discrepancy  could  only  be  resolved  by 
constructing  a variable  stress  drop/variable  rupture  veloc- 
ity model  for  the  faulting.  We  believe  the  same  conclusion 
is  warranged  here. 

We  construct  a second  model  for  the  Pocatello  Valley 
earthquake.  Model  II.  Our  hypothesis  is  that  the  earthquake 
must  have  included  regions  of  faulting  with  relatively  low 
stress  drop  and  rupture  velocity  that  are  important  for  the 
long  period  radiation,  but  are  only  weak  radiators  of  the 
high  frequency  radiation  seen  on  the  short-period  recordings. 
Many  combinations  of  fault  parameters  were  tried.  The  "best" 
of  these  will  be  described.  The  solution  is,  of  course,  far 
from  unique. 

The  faulting  for  Model  II  includes  two  regions.  In 
the  first  region  the  source  parameters  are  constant  and  not 
unlike  the  values  for  Model  I.  The  remainder  of  the  faulting 
is  characterized  by  a rupture  velocity  that  decays  to  zero, 
linearly  with  rupture  time.  In  this  region  the  stress  drop 
decays  to  zero,  but  with  a cubic  dependence  on  rupture  veloc- 
ity. This  rapid  decay  of  stress  drop  was  chosen  to  compen- 
sate for  the  cubic  dependence  of  amplitude  on  source  dimension 
that  is  peculiar  to  the  Archambeau/Mins ter  model.  A less 
rapid  decay  of  stress  drop  generated  too  much  long  period 
energy  in  the  short  period  synthetics.  In  summary,  we  have 
the  following: 

MODEL  II 

Fault  Length,  L:  7 km 

Rupture  Velocity,  V : 3 km/sec,  0 < L < 4 

^ 1/2 

(21-3L)  km/sec,  4 < L < 7 

Stress  Drop,  Aa : 46  bars,  0 £ L £ 4 

1.7  (21-3L)^'''^  bars,  4 < L < 7 

(unilateral  fault  propagation  toward  the  surface) 


Depth/  H;  8.7  km 

Orientation:  Strike  N45°,  Dip  39 °W 

Slip  Up  SS’N 

The  dependence  of  cind  Aa  on  position  along  the  fault  is 
plotted  in  Figure  9. 

As  far  as  the  long  period  records  are  concerned. 

Models  I and  II  are  essentially  the  same.  The  short  period 
synthetics  for  both  models  appear  in  Figure  8.  The  Model 
II  records  are  clearly  longer  period  and  are  missing  some 
of  the  high  frequency  detail  exhibited  in  the  Model  I records. 
As  far  as  the  major  phases  are  concerned,  the  P and  pP,  the 
Model  II  records  seem  preferable  at  some  stations  while  at 
others  no  advantage  is  apparent  with  either  model.  Another 
attractive  feature  of  Model  II  is  that  its  dimension  is  more 
consistent  with  the  size  of  the  region  defined  by  the  after- 
shock locations  (Figures  2 through  4) . Perhaps  a better 
model  could  be  found  somewhere  between  the  two;  they  seem 
to  span  a range  of  desirable  models. 

The  amplitude  data  for  Model  II  are  tabulated  in  Table 
6 in  the  same  format  as  the  Model  I data  in  Table  5.  We  see 
that  the  amplitude  ratios  are  2 to  3 times  larger  for  Model 
II  than  they  were  for  Model  I.  While  this  model  has  the  cor- 
rect long  period  level,  the  synthetic  short  period  records 
are  still,  on  the  average,  about  a factor  of  two,  larger 
than  the  observations.  It  may  be  that  this  is  as  close  as 
it  is  reasonable  to  get,  given  the  many  errors  inherent  in 
the  procedure  and  the  relatively  poor  quality  of  the  data. 

Attempts  to  improve  the  agreement  must  focus  on  con- 
structing a source  that  has  an  even  larger  region  that  con- 
tributes to  the  long  period  energy  without  substantially 
altering  the  short  periods.  This  is  difficult  because  the 
long  period  level  is  constrained  by  data  at  periods  that 
are  not  all  that  long;  five  or  six  seconds.  Also,  the  short 
period  waveforms  for  Model  II  are  already  too  long  period. 
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Another  interesting  and  potentially  important  aspect 
of  these  data  is  that  they  show  evidence  of  a small  precursor 
event  a second  or  so  before  the  main  event.  On  the  original 
Develocorder  records  for  each  of  the  five  stations  studied 
there  is  a small  arrival  just  before  the  main  P arrival.  Un- 
fortunately, our  digitization  is  too  poor  to  see  this  very 
clearly  in  the  Figure  8 plots,  though  it  can  be  discerned  on 
the  RES,  DAG  and  PTO  playouts.  The  sharp  break  in  the  wave- 
form between  the  major  P and  pP  arrivals  is  consistent  in 
amplitude  and  time  with  a pP  from  such  a precursor  event  at 
roughly  the  same  depth  and  orientation  as  the  main  event. 

This  event  appears  to  be  about  20  percent  as  large  as  the 
main  event  at  1 to  2 Hz  and  to  be  at  approximately  the  same 
depth  and  orientation.  The  time  lag  of  from  1 to  1.5  seconds 
could  be  a result  of  both  spatial  and  temporal  separation. 

2.6  CONCLUSIONS 

We  first  summarize  the  way  we  have  developed  our  model 
for  the  1975  Pocatello  Valley  earthquake.  We  have  a rather 
well-constrained  focal  mechanism  solution  for  the  event. 
Choice  between  the  two  original  solutions  was  made  on  the 
basis  of  the  depth  distribution  of  aftershock  locations  de- 
termined by  Arabasz,  e_t  al.  (1975b).  Thus,  we  begem  our  syn- 
thetic seismogram  studies  with  considerable  confidence  in  our 
knowledge  of  the  fault  orientation,  strike,  dip  and  slip. 

The  long  period  P waves  are  sensitive  to  the  fault 
orientation,  focal  depth  and  moment.  The  orientation  and 
depth  control  the  waveform  while  the  moment  controls  the 
amplitude.  Our  model  fits  the  first  5 to  8 seconds  of  the 
long  period  observations  (Figure  7) . The  moment  we  infer  is 
essentially  the  same  as  that  obtained  by  independent  investi- 
gation using  other  data.  Later  phases  on  the  long  period 
records  are  evidence  for  the  complexity  of  the  faulting 
after  the  first  several  seconds. 


The  main  phases  on  the  short  period  observations  are 
the  direct  P from  fault  initiation  and  the  associated  pP. 

The  sP  phase  is  small  for  this  fault  orientation.  These 
data  are  not  of  ^ood  quality  and  it  is  difficult  to  make 
detailed  comparisons  between  observed  and  synthetic  records. 
The  main  diagnostic  features  are  the  waveform  and  frequency 
content  of  the  P and  pP  phases,  their  spacing  and  their 
relative  amplitudes. 

We  construct  Model  I,  a single  rupture  velocity, 
single  stress  drop  model  for  the  earthquake  that  gives  a 
good  qualitative  fit  to  the  observations  (Figure  8) . One 
discrepancy  is  that  the  synthetics  appear  to  be  a bit  too 
short  period.  However,  if  this  source  is  to  have  the  right 
moment  (inferred  from  the  long  period  data) , the  synthetic 
short  period  records  are  much  too  large — a factor  of  3 to  5. 
We  see  no  way  to  resolve  this  discrepancy  except  to  suppose 
that  the  event  includes  regions  that  contribute  to  the  long 
period  level  without  substantially  adding  to  the  short  period 
radiation.  That  is,  we  need  a variable  stress  drop,  variable 
rupture  velocity,  model. 

We  construct  Model  II  which  begins  like  Model  I,  then 
moves  into  a region  where  the  stress  drop  and  rupture  veloc- 
ity decay  to  zero.  This  model  substantially  reduces  the 
discrepancy  between  the  source  amplitudes  inferred  from  the 
long  and  short  period  data.  The  main  problem  with  Model  II 
is  that  in  decreasing  the  short  period  amplitude,  we  have 
made  the  synthetic  seismograms  too  long  period.  That  is,  the 
variable  rupture  velocity,  variable  stress  drop  effects  are 
still  not  strong  enough  to  give  an  ideal  fit  to  the  data. 
However,  given  the  rather  large  errors  inherent  in  the  proce- 
dure, especially  in  comparing  amplitudes  at  periods  where  the 
frequency-dependent  instrument  response  is  rapidly  changing. 
Model  II  can  be  said  to  give  a reasonably  good  fit  to  both 
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the  long  and  short  period  data.  The  fit  could,  probably,  be 
improved  but  an  improved  model  must  have  roughly  the  same 
variable  stress  drop,  variable  rupture  velocity  character  as 
our  Model  II. 

The  stress  drop  at  initiation  of  Model  II  is  45  bars. 
This  compares  to  the  value  of  206  bars  found  with  the  same 
methods  for  the  1971  San  Fernando  earthquake  (Bache  and 
Barker,  1978) . In  our  discussion  of  that  event  we  point  out 
that  the  Archambeau/Minster  source  model  consistently  under- 
estimates the  stress  drop  by  a factor  of  about  3.6.  Thus,  we 
prefer  to  estimate  the  actual  stress  drop  for  the  initial 
faulting  to  be  about  162  bars.  In  any  case,  this  is  a factor 
of  4-1/2  smaller  than  the  San  Fernando  stress  drop  estimate. 

The  meaning  of  this  stress  drop  difference  is  diffi- 
cult to  assess  and  requires  further  thought.  It  is  fair  to 
say  that  all  stress  drop  estimates  are  "average"  estimates 
in  some  sense.  It  may  be  that  the  Pocatello  Valley  event 
includes  smaller  regions  with  substaintially  higher  stress 
drop  that  are  not  resolvable  with  the  data  at  hand.  Still, 
we  feel  confident  that  the  difference  is  qualitatively  cor- 
rect and  that  the  Pocatello  Valley  earthquake  was  a lower 
stress  drop  event. 

In  the  next  section  we  discuss  an  entirely  different 
method  for  computing  the  earthquake  source.  This  is  by  three- 
dimensional  finite  difference  modeling  of  stress  relaxation 
on  a fault  plane. 
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III.  THREE-DIMENSIONAL  FINITE  DIFFERENCE  EARTHQUAKE 

MODEL 


3.1  INTRODUCTION 

A three-dimensional  finite  difference  model  CTRES)  for 
earthquake  faulting  (Cherry,  1977)  has  been  made  operational 
on  the  parallel  processing  ILLIAC  IV  computer.  As  part  of 
this  contract,  ILLIAC  IV  computer  time  was  provided  to  model 
earthqucUce  faulting  with  particular  attention  to  the  Pocatello 
Valley  event.  In  this  section  we  present  the  results  of  our 
earthquake  modeling  with  the  ILLIAC  program. 

While  the  TRES  algorithm  is  quite  flexible  in  the  al- 
lowed specification  of  geometry  and  material  behavior,  the 
version  currently  operational  on  the  ILLIAC  is  quite  restricted. 
(The  allowed  fault  plane  geometry  is  shown  in  Figure  11.)  The 
faulting  must  be  bilateral;  in  fact,  it  must  nucleate  from  a 
point  and  propagate  with  circular  symmetry  until  reaching  the 
edges  of  a rectangular  fault  plane.  There  are  no  material 
boundaries  other  than  the  fault  plane;  that  is,  the  calcula- 
tions are  done  in  a whole  space. 

The  material  behavior  is  linearly  elastic  except  in 
the  vicinity  of  the  fault  plane  where  plastic  yielding  is 
permitted.  Including  this  rather  elementary  form  of  inelastic 
material  behavior  represents  an  important  first  step  toward 
developing  realistic  models  for  the  true  physics  of  earth- 
quake faulting. 

Two  full  calculations  were  carried  out  and  the  results 
will  be  described  in  this  section.  For  both  calculations  the 
material  properties  and  geometry  are  close  to  those  inferred 
for  the  Pocatello  Valley  earthquake  in  the  previous  section. 

The  difference  between  the  two  is  that  the  elastoplastic  be- 
havior near  the  fault  plane  was  allowed  in  one  case  and  not 
in  the  other. 


In  this  section  we  describe  the  finite  difference 
modeling  in  detail  and  present  the  results  as  they  apply  to 
earthqucukes  in  general.  In  later  sections  we  will  relate 
these  results  to  the  Pocatello  Valley  earthquake  and  discuss 
the  implications  for  western  U.S.  earthquakes  in  more  general 
terms . 

Section  3.2  describes  the  problem  formulation,  and 
Section  3.3  gives  the  fault  and  material  parameters  employed 
in  the  two  finite  difference  calculations.  In  Section  3.4, 
we  discuss  methods  for  continuing  the  numerical  solutions  to 
distances  beyond  the  calculational  mesh.  . In  Section  3.5,  we 
compare  our  numerical  solution  for  the  purely  elastic  case  to 
available  analytical  and  numerical  results  for  similar  prob- 
lems, in  order  to  verify  the  accuracy  of  the  numerical  method. 
Section  3.6  presents  the  results  for  the  near  source  stress 
and  velocity  fields  for  both  finite  difference  calculations. 
The  far-field  radiation  for  both  calculations  is  discussed  in 
Section  3.7. 

3.2  PROBLEM  FORMULATION 

We  treat  faulting  as  a propagating  stress  relaxation 
due  to  shear  failure  on  a planar  surface.  Ideally,  we  would 
like  to  specify  the  relevant  physical  properties  of  the  medium 
and  it.j  initial  conditions,  and  allow  a mathematical  model  of 
failure  to  determine  the  subsequent  evolution  of  the  fault 
plane.  However,  this  study  does  not  address  the  physical 
mechanism  of  failure.  Instead,  we  prescribe  the  propagation 
of  the  fault  surface.  Boundary  conditions  on  the  fault  sur- 
face are  governed  by  a simple  Coulomb  friction  law. 

We  will  specify:  (i)  the  initial  state  of  the  medium 
and  its  constitutive  properties,  (ii)  the  evolution  of  the 
fault  plane,  and  (iii)  the  boundary  conditions  to  be  satis- 
fied on  the  fault  plane.  The  mathematical  model  of  faulting 


is  then  solved  using  a three-dimensional  finite  difference 
method  (Cherry,  1977). 


3.2.1  Initial  Conditions  and  Constitutive  Properties 

For  time  t less  than  zero  we  assxime  that  an  equilib- 
rium state  of  stress  exists,  with  displacement  and  velocity 
everywhere  zero.  The  equilibrium  configuration  is  such  that 
the  prospective  fault  plane  experiences  a shear  stress 
and  a compressional  normal  stress 

Average  stress  changes  associated  with  faulting  are 
modest  — on  the  order  of  a few  hundred  bars.  Furthermore, 
faulting  represents  a relaxation  of  stress,  except  near  the 
fault  edges.  Linear  elasticity  is  thus  probably  am  adequate 
model  of  material  behavior  outside  the  immediate  zone  of 
faulting.  However,  large  strains  do  occur  at  the  edges 
of  the  fault  ahead  of  the  rupture.  Some  model  of  rock  strength 
must  be  incorporated  to  prevent  large  stress  concentrations 
from  accruing  immediately  ahead  of  the  fault. 

For  this  study,  a simple  model  of  plastic  yield  was 
utilized  in  the  region  adjacent  to  the  fault.  The  plastic 
flow  model  is  as  described  in  Cherry,  et  al.  (1976).  If  a' 
is  the  deviatoric  stress  tensor,  calculated  assuming  that  the 
strain  rate  is  elastic,  then  the  actual  stress  deviator  a' 
is  given  by 


Y" 

^2  > T 


for  >^2  - 3“  ' 


where  is  the  second  deviatoric  invariant  defined  by 
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and  Y is  a material  constant  representing  yield  strength.  The 
stress  adjustment  described  by  Equation  (3)  is  permitted  only 
within  a specified  zone  near  the  fault.  Elsewhere,  linear 
elasticity  is  presumed  to  hold. 

3.2.2  Evolution  of  the  Fault  Pletne 

For  this  study,  we  specify  the  geometry  of  the  fault 
surface  as  a function  of  time,  rather  than  deriving  the  fault 
surface  from  the  model  as  a consequence  of  a failure  mechanism. 
This  is  viewed  as  a first  step  toward  more  fully  deterministic 
calculations . 

The  rupture  surface  is  assumed  to  be  plauiar,  to  nucle- 
ate from  a point,  and  to  expand  symmetrically  at  a constant 
prescribed  rupture  velocity  (V  ) , until  it  reaches  prescribed 
boundaries.  SCt)  denotes  the  fault  surface.  If  r and  9 are 
polar  coordinates  in  the  fault  plane,  with  r being  the  hypo- 
central  distance,  and  if  the  edge  of  E (<»)  is  defined  by  the 
curve 


r = 8(0)  , 

then  1 t)  consists  of  all  points  r,  9 such  that 
r < min  [Vj^t,B(9)]  H(t)  . 

Two  significant  features  of  this  model  of  rupture  are 
(i)  the  rupture  velocity  accelerates  instantaneously  to  its 
final  velocity,  and  (ii)  the  rupture  advance  terminates 
abruptly,  i.e.,  the  rupture  velocity  decelerates  instanta- 
neously to  zero  at  some  prescribed  boundary.  The  former  as- 
sumption may  be  quite  reasonable.  There  is  both  theoretical 
(Cherry,  et  al. , 1976;  Das  and  Aki , 1977)  and  experimental 


(Archuleta  and  Brune,  1975)  evidence  that  ruptures  accelerate 
rapidly  to  their  terminal  velocity,  with  the  rupture  velocity 
essentially  constant  over  most  of  the  fault.  The  latter  as- 
sumption, that  stopping  of  the  rupture  is  abrupt,  is  probably 
an  unrealistic  model  for  many  earthquakes.  For  example,  Bache 
and  Barker  (1977)  infer  a fault  model  for  the  San  Fernando 
earthquake  in  which  the  rupture  velocity  decelerates  smoothly 
to  zero.  On  the  other  hand,  if  the  rupture  termination  is 
controlled  by  a fracture  energy  barrier,  abrupt  stopping  may 
be  appropriate  (Husseini,  ^ al. , 1975). 

3.2.3  Boundary  Conditions  on  the  Fault 

On  E (t)  we  permit  a tangential  displacement  discon- 
tinuity £(x,t) , and  require  continuity  of  the  normal  displace- 
ment: 


£(x,t)  = lim  [u(x  + en)  - u(x  - en)  ] , x on  E, 

e-^0  (4) 

A 

n • s = 0 , 

with  £ > 0.  The  n is  a unit  vector  normal  to  E and  u is  the 
displacement  vector.  We  also  require  continuity  of  traction. 

The  tangential  traction  on  E is  assumed  to  obey  a 
simple  Coulomb  friction  law.  Let  denote  the  tangential 
traction  vector,  defined  in  terms  of  the  stress  tensor  a by 

jr  s - n * a • (^  - M)  , (5) 

where  ^ is  the  indentity  tensor.  Then  define  to  be  a slid- 
ing frictional  traction  which  opposes  the  slip  velocity  s^  and 
is  proportional  to  the  normal  stress: 


> 0, 


where  is  - u n • a • h/  the  product  of  the  normal  stress 
and  the  coefficient  of  dynamic  friction  (a ^ is  presumed  to  be 
positive).  Finally,  define  as  the  tangential  traction  on 
E which  would  be  required  to  enforce  continuity  of  velocity. 

We  wish  to  set  the  tangential  traction  equal  to  T_f 

whenever  the  slip  velosity  ^ is  non-zero.  At  the  instant  when 

• 

s becomes  zero  at  a given  point,  equals  x , by  definition. 

If  further  slip  is  permitted  at  that  point,  t_  reverse  to 

However,  there  is  a condition  under  which  further  slip 
is  inconsistent  with  the  equation  of  motion.  To  see  this,  we 
write  the  equation  of  motion  (in  the  absence  of  body  forces) 


pu  = V 


Letting  the  fault  be  the  plane  z = 0 , we  consider  a small 
sphere  surroionding  the  point  x on  the  fault.  Figure  10  il- 
lustrates this.  and  S denote  the  surface  of  the  hemi- 

sphere in  the  +z  and  -z  halfspaces,  respectively.  Taking 
the  scalar  product  of  Equation  (7) , first  with  ^ sgn  (z) , 
then  with  y sgn  (z) , integrating  over  the  volume  of  the 
sphere,  and  applying  the  divergence  theorem,  we  get  the  pair 
of  equations : 

::  - ~ 

""x  = -FT  (""x  ^ -x)  ' 


Au  = (Af  + T ) . 

y M y y 
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Aj,  is  the  area  of  the  fault  plane  encompassed  by  the  sphere; 

M is  the  mass  of  each  hemisphere;  Au  and  Au  are  the  x amd  y 

X y 

components,  respectively,  of  the  average  acceleration  in  the 
+z  hemisphere  minus  the  average  acceleration  in  the  -z  hemi- 
sphere; Af  and  Af  are  the  x and  y components,  respectively, 

^ y ^ 

of  the  average  traction  on  the  exterior  of  S minus  the  aver- 
age traction  on  the  exterior  of  S ; and  x^  are  the  compo- 
nents, respectively,  of  x,  averaged  over  A^. . As  the  radius 
of  the  sphere  is  reduced,  Au^  and  Au^  approach  s^  and  s^ , 
respectively;  and  x^  and  x^  approach 


■ ""f  / 


/•2  ^ *2  \ 
Is  + s 1 
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respectively,  if  slippage  is  occurring.  Thus,  Equation  (9) 
takes  the  form 


2Aj.  ^ s 

= — “x  - “£  jr-|5-p[72-  ■ 

X y 


- M - “f  JirrtTW^  • 


If  slip  is  to  recommence  after  going  to  zero  at  t = t^,  s_  (t^) 

will  have  the  same  direction  as  £ (t"*”)  , so  the  equations  at 

+ ° 
t can  be  written 
o 


lAf 

s..  = -rr-  Af..  - a 


M - '^f  (■■2  ^ ^2072 


(10a) 


(10b) 


This  pair  of  non-linear  equations  in  two  unknowns  has  the  solu 
tion  (Day,  1977,  Appendix  IV) 


s..  = 


M 


Af..  - a. 


sgn  (s^) 


(11a) 


s = 


2^ 

M 


Af  - a 


sgn  (^) 
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(11b) 


This  solution  would  always  exist  if  the  coefficient  of  the 
signum  function,  a^,  were  negative.  Note,  however,  that  this 
coefficient  is  necessarily  positive,  by  virtue  of  the  fact 
that  friction  always  opposes  the  slip  velocity.  As  a result, 
there  exists  a condition  under  which  no  solution  exists  to 
the  system  given  in  Equation  (10) . The  condition  for  no 
solution  is 


> Af2\l/2 


< 0 


(12) 


This  is  just  the  condition  that  exceeds  |t_^|  , in  which 
case  we  see  that  further  slip  is  precluded.  This  leads  to 
the  following  statement  of  the  boundary  conditions  on  E govern 
ing  T^: 


If,,  if 


> 


T = 


if 


llcl  1 


(13) 


Thus,  healing  of  the  fault  is  not  arbitrarily  enforced,  but  re 
suits  when  the  equations  of  motion  will  not  admit  a solution 
with  a velocity  discontinuity. 
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3.3  FAULT  AND  MATERIAL  PARAMETERS 


We  performed  two  three-dimensional  finite  difference 
calculations  for  this  study.  The  two  calculations  differed 
only  in  the  yield  strength  Y assigned  the  material.  Both  cal- 
culations were  for  a square  fault  plane  in  a uniform  whole 
space,  with  rupture  initiated  at  the  center  of  the  square 
fault  (Figure  11) . The  following  parameters  were  employed 
for  both  calculations: 


P wave  velocity 
S wave  velocity 
Density 

Rupture  velocity 
Tectonic  shear  stress 
Frictional  stress 
Fault  dimensions 


a =5.93  km/sec 

B =3.42  km/sec 

p =2.74  gm/cm 

V„  =3.08  km/sec 

R 

=1  kbar 
=0.82  kbars 
2ax2a  = 3 km  x 3 km 


For  the  first  finite  difference  calculation,  the 
"elastic"  model,  the  yield  strength  Y was  set  to  infinity, 
so  the  constitutive  model  was  linearly  elastic.  For  the 
second  finite  difference  calculation,  the  "plastic"  model, 

Y was  set  to  /3a^,  so  that  the  fault  zone  was  initially 
stressed  to  the  failure  surface.  With  this  choice  of  Y, 
plastic  flow  ensues  when  the  deviatoric  invariant  increases 
above  its  initial  value.  This  dissipates  any  dynamic  shear 
stress  concentration  ahead  of  the  crack  tip.  Plastic  yield 
was  permitted  only  within  0.2  km  of  the  fault  plane;  else- 
where linear  elasticity  was  employed. 

For  both  the  elastic  and  the  plastic  fault  problems, 
the  medium  was  represented  by  cubic  finite  difference  zones 
0.1  km  on  a side.  The  nvunerical  grid  was  large  enough  that 
no  reflection  from  the  exterior  grid  boundary  returned  to 
the  fault  zone  during  the  calculation. 
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The  formulation  of  Section  3.2  was  simplified  somewhat  for 
these  finite  difference  calculations  by  prohibiting  fault 
slip  in  the  direction  normal  to  the  prestress  direction. 

This  restriction  should  not  have  a significant  effect  on  the 
solution.  It  is  known,  for  example,  that  for  the  self-similar, 
expanding  elliptical  crack,  the  slip  is  everywhere  parallel 
to  the  prestress.  Furthermore,  Madariaga  (1976)  has  shown 
that  the  perpendicular  component  of  slip  is  very  small  for 
the  case  of  a finite  circular  crack. 

3.4  CONTINUATION  OF  THE  ELASTIC  FIELD 

The  computing  cost  for  three-dimensional  finite  dif- 
ference calculations  increases  rapidly  with  the  distance  over 
which  waves  must  be  propagated  through  the  numerical  grid. 

Since  we  are  interested  in  simulating  earthquake  ground  motion 
at  a variety  of  distances,  it  is  important  to  be  able  to  link 
the  numerically  computed  wave  fields  to  analytic  wave  propaga- 
tion methods.  A very  general  method  for  doing  so  is  provided 
by  an  elastodynamic  integral  representation  given  by,  for 
example,  Burridge  and  Knopoff  (1964).  The  approach  is  approp- 
riate in  both  the  near-field  and  far-field.  Two  elements  are 
necessary  for  this  analytic  continuation:  (i)  Time  histories 
of  the  stress  vector  x_  (^,t)  and  the  displacement  vector 
u (^,tj  are  required  on  a surface  Z surrounding  the  source 
(and  enclosing  any  inelastic  region) . These  are  available  from 
the  finite  difference  source  calculation.  (ii)  The  Green's 
(tensor-valued)  function  g spatial 

derivatives  are  needed  for  the  elastic  medium  representing  the 
propagation  path  to  the  observation  point  of  interest  x.  For 
a horizontally  layered  earth  structure,  these  can  be  readily 
computed  by  wave  propagation  techniques . 

When  (i)  and  (ii)  are  evaluated,  the  solution  for  mo- 
tion at  a given  site  is  obtained  via  spatial  integrals  over 
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the  surface  surrounding  the  source.  The  general  form  of  the 
integral  representation  (assuming  isotropy)  is: 


00 

u (x,t)  ~ f f 


dS  [g 


- g V (x,x  ; t-t  ) : M (x  »t  )] 
* — — — o o = — o o 


(14) 


where 


= ^In  • u (x  ,t  ) + y [nu  (x  ,t  ) 
**  — o o — o o — o o 


+ u (x  ,t  ) n]  , 

- -o  o - 


X and  t are  receiver  position  and  time, 

X and  y are  elastic  constants, 
n IS  the  normal  to  E,  directed  into  the  source  volume,  * i 

u is  the  displacement  vector,  .J 

I 

< 

T_  is  the  stress  vector  on  E , 1 

g is  the  Green's  tensor  solution. 


^ is  the  identity  tensor. 

In  the  special  case  in  which  inelasticity  near  the  fault  is 
negligible,  the  surface  can  be  reduced  to  the  fault  plane. 

In  this  case,  the  continuation  method  reduces  to  the  so- 
called  dislocation  method  (see,  for  example,  Has)cell,  1964; 
Aki,  1968)  — only  the  displacement  discontinuity  at  the 
fault  need  be  monitored  in  order  to  propagate  the  source 
disturbance. 

Two  special  cases  of  the  representation  are  particu- 
larly useful  for  isolating  the  far-field  radiation.  The 
first  method  applies  to  any  spherical  source  region  embedded 
in  a uniform  elastic  whole  space.  The  second  applies  to  a 
planar  failure  surface  embedded  in  a uniform  whole  space. 
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We  designate  these  the  "multipolar  expansion"  method  and  the 
"angular  spectrum"  method,  respectively.  We  will  apply  these 
methods  in  subsequent  sections  to  obtain  the  radiated  seismic 
field  from  the  finite  difference  solutions. 


3. 4. 1 Multipolar  Expansion  Method 

When  the  source  region  is  embedded  in  a homogeneous 
whole  space,  we  can  simplify  the  representation  of  the  ex- 
terior seismic  field.  This  is  accomplished  by  choosing  the 
integration  surface  to  be  a sphere,  writing  the  displacements 
and  Green's  functions  in  terms  of  scalar  and  vector  potentials 
(divergence  and  curl)  and  expanding  the  scalar  potentials  as 
spherical  harmonic  series.  The  approach  applies  equally  well 
to  the  elastic  and  elastoplastic  cases  in  this  study,  pro- 
vided the  integration  sphere  surrounds  the  region  of  non- 
linearity. Bache  and  Harkrider  (1976)  and  Bache,  et  al . 

(1975)  give  details  of  the  formulation.  The  following  is  a 
brief  outline  of  the  method. 

The  Fourier  transfoirmed  equations  of  motion  in  a homo- 
geneous, isotropic,  linear  elastic  medium  may  be  written 


u = - 


Yx 


(4) 


Y X X 


(15) 


a; 


where  u is  displacement  and  k and  kp,  are  the  compressional 


and  shear  wave  numbers.  The  Cartesian  potentials  x 
X are  defined  by 


(4) 


and 


x^^^  = Y • u , 


(16) 


X = - V X u 


The  potentials  can  be  shown  to  satisfy  wave  equations  and 
therefore  can  be  expanded  in  spherical  eigenfunctions  as 
follows : 


(x,aj) 


= E ■'f’  E [' 


(to)  cos  mi)) 


1=0 


m=0 


+ (u)  sin  (cose),  j = 1,2, 3, 4 (17) 


where  k.  h = to/a  and  k.  H k„  = to/6  for  i = 1,2,3.  The 

( 2 ) ^ ^ 1 p 

h'  are  spherical  Hankel  functions  of  the  second  kind  and  the 
are  associated  Legendre  functions,  x is  the  position  vec- 
tor, and  r,0,i))  are  spherical  coordinates  such  that  x is  given 
in  terms  of  Cartesian  unit  vectors  x,  £,  £ by 


X = X r sine  cose))  + y r sine  sini))  + £ r cose  . 


The  multipole  coefficients  Af^^(to),  Bf^^(to),  j = 1,..., 

x/in  Join 

4,  are  determined  from  the  observed  divergence  and  curl  of  the 
displacement,  obtained  from  the  finite  difference  solution. 
Using  the  orthogonality  of  the  spherical  harmonics,  the  multi- 
pole coefficients  are  given  by  the  expression 


^(j) 

ilm 
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(j) 

ilm 


ilm 

h/’  H'c.R) 


TT 


1 cos  mij)  j 

(x,a))  P^(cose)  I ( 

sin0dedi)) 

( sin  mi)); 

(18) 

where 


C 


ilm 


(2il-H)  (il-m)  i 
2tt  ( il+m) 


m / 0 


C,Q  = (2^+1)/4tt, 

and  R is  the  radius  of  the  spherical  surface  on  which  the 
integral  is  to  be  computed. 


The  multipole  coefficients  Aj|^^  (to)  , B^^Nw)  , j = 1,..., 
4 specify  the  displacement  field  (through  Equations  15  and  17) 
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everywhere  exterior  to  the  (spherical)  source  region.  This 
is  a convenient  form  for  studying  whole  space  radiation 
properties.  Separation  into  near-  and  far-field  components 
is  simple,  since  the  r dependence  is  isolated  in  the  spheri- 
cal Hankel  functions.  Also,  radiation  patterns  are  easily 
exhibited,  since  the  angular  dependence  is  isolated  in  the 
spherical  harmonics. 

One  significant  drawback  to  this  approach  is  the 
necessity  to  compute  the  full  transient  finite  difference 
solution  on  the  spherical  surface.  This  requires  considercibly 
more  computation  time  than  is  required  to  complete  the  tran- 
sient solution  on  the  fault  plane.  A second  difficulty  is 
that  there  is  some  uncertainty  in  the  number  of  terms  which 
must  be  retained  in  the  expansion  in  order  to  obtain  conver- 
gence at  a given  frequency.  Finally,  the  expansion  is  poorly 
convergent  when  the  distance  from  the  origin  to  the  observa- 
tion point  is  comparable  to  the  source  dimension. 

3.4.2  Angular  Spectrum  Method 

When  non-elastic  behavior  can  be  neglected,  the  inte- 
gration surface  can  be  collapsed  to  the  fault  plane.  For 
this  case,  Kostrov  (1968)  and  Dahlen  (1974)  review  the  result- 
ing simplification  of  Equation  (14)  for  finding  the  far-field, 
whole  space  radiation.  The  first  term  in  Equation  (14) 
vanishes,  due  to  the  continuity  of  traction  at  the  fault 
plane.  Further  simplification  results  from  employing  the  far- 
field  approximation  to  the  whole  space  Green's  tensor,  the 
point-source  approximation  (source  dimension  much  less  than 
source-receiver  distance) , and  from  assuming  a planar  fault 
surface  with  uni-directional  slip  tangential  to  the  fault 
plane.  This  leads  to  Dahlen 's  (1974)  Equations  (13)  to  (16) 
for  the  Fourier  transformed  P and  S wave  displacements  u'^  and 

Q ■“ 

u^ , which  can  be  written  as 


-i  iiir/a 


u“  (r,  6 ,({)  ,(jj)  = " 


R (9,(j))  n (0,4), (d) 


47t  a r 


6 e"^  - 

u (r,e,4>,a))  = — gj.-  [6^3  0 


+ (Q/'J')]  (6/<{)/(jj), 

where  R , R „ , and  R . are  the  double  couple  radiation  pat- 
tern  components.  They  are  given  in  terms  of  n,  e and  r,  the 
unit  vectors  in  the  directions  of  the  fault  normal,  slip 
direction,  and  source- receiver  direction,  respectively,  by 

/S  ^ /S  /V  A. 

R = r r:  (n  e + e n), 

p - - _ - _ _ , 

(20: 

R g = § r:  (n  e + e n). 


R , = 4)  r:  (n  e + e n) 
s4)  - ~ ~ ~ - 


For  the  problem  geometry  of  Figure  11,  these  factors  are 


R = sin  29  sin4). 
P 


R „ = cos  29  sin4», 


R.  . = cos9  cos4)  . 
s4i 


The  fJp  and  are  obtained  directly  from  a three-dimensional 

Fourier  transfoirm  (the  "angular  spectrum"  (Bracewell,  1965)) 

of  the  slip  velocity  function.  TaJcing  x^,  y^  to  be  Cartesian 

coordinates  describing  the  fault  plane,  and  denoting  the 

velocity  discontinuity  across  the  fault  by  e s (x^,y^,t), 

n and  are  given  by 
p s 


L.. 
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(9  »<))  /OJ) 
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^ sin6  cos4> 


^ sin6  sin(() 
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Jig  (0  » /CO) 
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(22a) 


j sin9  cos<j> 


Ky  = j sin9  sin<p 
(22b) 


Equations  (22a)  cuid  (22b)  contain  all  the  dependence 
of  the  radiation  on  the  finite  difference  solution.  Thus,  it 
is  only  necessary  to  store  the  computed  slip  velocity  at  all 
nodes  on  the  fault  pleme,  and  it  is  only  necessary  to  com- 
pute the  finite  difference  solution  until  the  fault  has  every- 
where healed.  Given  this  slip  velocity  function,  Equations 
(22a)  and  (22b)  require  a single  three-dimensional  Fourier 
transform  be  performed.  Then  (or  n^)  as  a function  of  co , 
for  any  given  9,4>,  is  obtained  by  assembling  the  Fourier 

transform  values  along  a line  through  K , K , co  space. 

X y 

Equation  (19)  then  gives  the  ? and  S wave  spectra. 

Since  s is  known  from  the  finite  difference  solution 
at  discrete  points,  it  is  natural  to  employ  an  FFT  algorithm 
to  execute  the  three-dimensional  transform.  The  application 
is  entirely  straightforward  and  efficient.  The  only  note 
of  caution  concerns  interpolation  of  the  discrete  transform. 

For  a given  ou,  determining  0.  requires  interpolating  the  K^, 

Ky.  dependence  of  the  discrete  transform.  The  finite  difference 
solution  provides  the  series  s.  . .,  i » j = 1,...J, 

k = 1,...  K,  where  I and  J are  the  number  of  nodes  in  the  x 
and  y directions,  respectively,  at  which  slip  is  non-zero, 
and  K is  the  number  of  time  steps  in  the  solution.  Since  the 
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function  s is  zero  outside  this  interval,  the  appropriate 
K , K interpolation  of  the  discrete  transform  is  via  convolu- 

^ y 

tion  of  the  transform  with  a two-dimensional  sine  function. 

This  is,  of  course,  best  accomplished  by  extending  the  series 
s by  adding  zeros  in  the  i and  j directions.  In  practice  it 
appears  sufficient  to  add  enough  zeros  to  double  the  series 
length  in  these  directions.  For  the  applications  presented 
in  subsequent  sections,  the  series  were  padded  to  four  times 
their  original  length  in  both  directions.  Further  interpo- 
lation in  Kx,  Ky  can  be  done  linearly  without  distortion  of 
the  solution. 

3.5  VERIFICATION  OF  THE  METHOD  AND  COMPARISON  TO  PREVIOUS 

WORK 

The  numerical  procedures  involved  in  solving  the  non- 
linear boundary  value  problem  (Sections  3.2  and  3.3)  and 
analytically  continuing  the  resultant  seismic  field  (Section 
3.4)  are  sufficiently  complex  that  a careful  evaluation  of 
their  accuracy  is  warranted.  The  ideal  verification  approach 
would  be  to  treat  a simple  case  for  which  an  exact  analytic 
solution  is  available,  then  compare  the  analytic  solution  to 
the  numerical  solution.  Unfortunately,  no  analytic  solution 
exists  for  a three-dimensional  problem  in  which  a stress  drop 
is  specified  on  a propagating  fault  surface  which  ultimately 
stops  growing  and  heals.  In  this  section  we  compare  our  numeri- 
cal solution  for  the  elastic  fault  problem  to  a variety  of 
available  analytic  and  numerical  results  for  similar  problems. 
Then  we  check  the  internal  consistency  of  the  method  by  com- 
paring the  far-field  pulses  obtained  from  Fourier  analysis  of 
the  slip  function  with  those  obtained  from  the  spherical 
harmonic  expansion  of  the  outgoing  elastic  field. 

3. 5. 1 Kostrov's  Analytical  Solution 

Kostrov  (1964)  obtained  cin  analytical  solution  for  the 
problem  of  a circular  crack  which  nucleates  at  a point  in  a 
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homogeneous,  unbounded  elastic  medium  and  expands  at  a con- 
stant rupture  velocity  without  stopping.  His  solution  for 
the  displacement  discontinuity  across  the  crack  is 


where  C is  a constant  which  depends  on  rupture  velocity  (see 
Dahlen,  1974) , is  che  excess  of  the  prestress  over  the 
kinetic  frictional  stress,  u is  the  shear  modulus,  3 is  the 
shear  wave  speed,  r is  distance  from  the  hypocenter,  and  V 

R 

is  the  rupture  propagation  velocity.  This  problem  corresponds 

exactly  to  the  condition  of  the  elastic  model  considered  in 

this  study,  until  time  sl/V^,  where  a is  the  half-length  of 

the  square  fault.  In  this  case  a = 1.5  km  and  V„  is  3.08 

R 

km/sec  (0.93),  for  which  C equals  approximately  0.81.  We 
shall  compare  the  initial  fault  slip,  obtained  numerically, 
to  Kostrov's  analytic  solution.  Once  the  rupture  front 
reaches  the  edge  of  the  square  fault  plcine  and  stops  growing, 
we  expect  the  numerical  solution  to  begin  to  deviate  sig- 
nificantly from  Kostrov's  solution. 

Figure  12  shows  the  slip  obtained  at  several  points  in 
the  fault  plane.  The  dashed  curves  are  the  finite  difference 
solution  and  the  solid  curves  are  the  analytic  solution.  The 
vertical  bars  indicate  the  arrival  times  of  edge  effects  due 
to  stopping  of  the  rupture  at  its  outer  boundary.  The  two 
solutions  display  the  anticipated  agreement  at  each  point 
prior  to  the  arrival  of  the  edge  effects.  The  small  devia- 
tion of  the  numerical  solution  from  the  analytic  solution  at 
early  time  results,  at  least  partially,  from  imprecise  weight- 
ing of  the  stress  drop  to  account  for  the  fractional  rupture 
of  a finite  difference  zone  by  the  circular  rupture  front. 
Archuleta  and  Frazier  (1978)  achieved  somewhat  better  agree- 
ment with  Kostrov's  solution  by  incorporating  fractional  rup- 
ture into  their  finite  element  scheme. 
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Figure  12.  Relative  displacement  on  the  fault  for  the  elastic 
case.  The  dashed  curves  are  Kostrov's  analytic 
solution;  the  solid  curves  are  the  finite  difference 
results.  x,y  coordinates  in  kilometers  are  given  in 
parenthesis.  Vertical  lines  indicate  the  arrival 
times  of  edge  effects  due  to  fault  finiteness. 


3.5.2  The  Stopping  Phase 


Once  edge  effects  occur,  we  have  to  rely  on  other  nu- 
merical solutions  for  comparison.  The  problem  of  a finite 
circular  fault  has  been  solved  by  Madariaga  (1976)  using  a 
finite  difference  method.  There  are  several  features  charac- 
terizing his  solution.  When  the  fault  stops  growing,  slip 
velocity  at  a given  point  on  the  fault  begins  to  deviate  from 
Kostrov's  solution  at  a time  corresponding  to  the  P travel 
time  to  that  point  from  the  nearest  edge  of  the  fault.  Thus, 
a stopping  phase  propagates  radially  inward  at  the  P wave 
speed.  The  subsequent  arrest  of  slip  also  initiates  at  the 
outer  edge  and  propagates  radially  inward.  The  velocity  of 
this  healing  "wave"  is  not  constant,  but  its  average  is  near 
the  shear  wave  velocity. 

Although  we  are  dealing  with  a square  fault  edge, 
similar  behavior  is  expected  since  rupture  growth  is  circu- 
lar and  our  final  fault  geometry  is  similar.  Figure  12  shows 
the  stopping  phase  for  points  on  the  square  fault.  The 
stopping  behavior  is  very  similar  to  that  of  the  circular 
fault,  in  spite  of  the  reduced  symmetry  in  the  problem  (the 
circular  fault  stops  everywhere  on  its  perimeter  simulta- 
neously, whereas  the  square  fault  does  not).  We  also  note 
that  t e displacement  curve  for  the  point  at  the  hypocenter 
differs  from  those  for  other  points  in  the  fault  plane  in 
that  it  has  two  distinct  breaks  in  slope  associated  with 
stopping.  This  feature  has  been  observed  in  several  numeri- 
cal studies  of  circular  shear  cracks  (Archuleta,  1976; 
Madariaga,  1976;  Day,  1977), 

3.5.3  Static  Configuration 

Neuber  (1937)  solved  the  static  problem  of  the  dis- 
placement of  an  elastic  whole  space  due  to  a uniform  shear 
stress  drop  on  a circular  crack  and  obtained  the  slip  distri- 
bution 


s (r) 


(24) 
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where  s is  the  relative  displacement  on  the  fault,  Aa  is  the 
stress  drop,  y the  shear  modulus,  b the  fault  radius  and  r 
distance  from  the  fault  center.  From  Equation  (24)  we  find 
that  the  average  slip  s is  2/3  the  slip  at  the  center  of  the 
fault.  However,  numerical  solutions  to  the  corresponding 
dynamic  problem  give  a maximum  and  average  static  slip  ap- 
proximately 25  percent  greater  than  Equation  (24) , for  the 
case  in  which  rupture  velocity  is  0.9  (interpreting  Aa  as  the 
difference  between  prestress  and  Icinetic  frictional  stress). 

For  example,  Madariaga  (1976)  found  a slip  function  20  per- 
cent in  excess  of  Equation  (24) , and  Archuleta  (1976)  and 
Day  (1977)  found  a 27  percent  overshoot.  The  failure  of  the 
dyncimic  solutions  to  approach  the  static  solution  is,  of 
course,  a consequence  of  the  non-linearity  of  the  boundary 
conditions  (i.e. , the  healing  of  the  fault) . 

From  the  square  fault,  we  obtained  a maximum  residual 
slip  of  129  cm,  which  occurred  at  the  center  of  the  fault. 

If  we  reinterpret  b in  Equation  (2  4)  as  /A/ir , where  A is  the 
fault  area,  then  129  cm  exceeds  the  prediction  of  Equation 
(24)  by  a factor  of  1.24,  in  good  agreement  with  numerical  re- 
sults for  the  circular  fault.  The  average  slip  for  the  square 
fault  was  79  cm,  which  is  0.61  times  the  value  of  the  center. 

The  seismic  moment  obtained  from  this  average  slip  was  2.28  x 
2 4 

10  dyne-cm.  This  value  is  14  percent  greater  than  the  pre- 
diction obtained  by  combining  the  static  crack  formula. 

Equation  (24)  , with  the  expression  for  seismic  moment,  = 
yAs. 

3.5.4  Computation  of  the  Radiated  Field 

As  a check  on  the  internal  consistency  of  the  analysis 
procedure,  we  compute  far-field  displacement  spectra  and  pulses 
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for  the  elastic  problem  by  two  different  methods.  Spectra  and 
pulses  obtained  from  the  angular  spectriom  of  the  slip  function 
are  compared  to  those  obtained  from  spherical  harmonic  expan- 
sions of  the  divergence  and  the  Cartesian  components  of  the 
curl,  evaluated  on  a sphere  of  2.7  km  radius.  Of  particular 
interest  is  the  convergence  rate  of  the  spherical  harmonic 
expansion,  since  we  will  rely  on  this  technique  alone  to  ob- 
tain the  radiation  for  the  plastic  crack  problem. 

Figure  13  compares  four  determinations  of  the  far- 
field  S wave  spectrum  and  time  pulses  for  the  elastic  fault 
model,  observed  at  9 = 0°.  The  dashed  curves  represent  the 
spectrum  and  pulse  obtained  by  Fourier  analysis  of  the  slip 
function.  The  other  three  sets  of  curves  were  obtained  by 
truncating  the  spherical  harmonic  series  at  Jl  = 4,  6,  and  8, 
respectively.  The  curves  have  been  normalized  by  the  zero- 
frequency  level  (4TTp8^r)  where  is  the  seismic  moment 
obtained  from  the  average  static  slip  s: 

= yAs  . (25) 

As  the  figure  indicates , the  spherical  harmonic  expansion  con- 
verges rapidly  at  low  frequency.  For  frequencies  above  1 to 
2 Hz,  however,  the  truncated  series  can  be  significantly  in 
error.  For  example,  at  3 Hz  the  spectrum  obtained  by  trun- 
cating the  series  at  £ = 6 is  in  error  (relative  to  the  spec- 
trum obtained  from  the  slip  function)  by  a factor  of  4. 
Including  terms  up  to  ?.  = 8 reduces  the  error  to  a factor  of 
1.5  (the  maximum  frequency  for  which  the  finite  difference 
solution  is  reliable  is  approximately  4 to  5 Hz) . 

Figure  14  compares  far-field  pulses  and  spectra  at  3 
azimuths;  the  dashed  curves  were  obtained  from  the  slip  func- 
tion, the  solid  curves  from  the  expansion  in  spherical 
harmonics  up  to  order  8.  For  computing  far-field  solutions 
in  the  range  0-3  Hz,  the  spherical  harmonic  expansion,  trun- 
cated at  = 8 appears  to  be  satisfactory. 
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Comparison  of  four  determinations  of  the  far-field 
S wave  spectrum  and  time  history,  for  the  elastic 
fault  model,  observed  at  6 = 0 . The  normalization 
is  described  in  the  text. 


Figure  14  . Far-field  time  histories  and  spectra  at  3 azimuths.  The  dashed  curves 
were  obtained  by  the  angular  spectrum  method,  the  solid  curves  by  the 
multipole  expansion  method  using  terms  up  to  order  £=  8.  The  normalization 
is  described  in  the  text. 
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3.6  NEAR-SOURCE  DEFORMATION 

In  the  last  section,  we  compared  the  slip  function  for 
the  elastic  fault  model  to  some  three-dimensional  analytical 
and  numerical  solutions  to  related  problems.  In  this  section, 
we  look  in  detail  at  the  near-source  stress  and  velocity 
fields  for  both  the  elastic  and  plastic  cases.  Our  main  ob- 
jective will  be  to  evaluate  the  effect  of  inelastic  material 
response  in  the  fault  zone. 

3.6.1  Stress  History  Near  the  Fault 

Figure  15  depicts  the  stress  histories  near  the  fault 
plane  for  the  two  models  (the  purely  elastic  model  and  the 
plastic  model) . The  shear  stress  in  the  direction  of  pre- 
stress plotted  versus  time  for  seven  points  along 

the  fault  plane  diagonal,  at  increasing  distance  from  the 
hypocenter.  The  stresses  shown  are  actually  evaluated  at  the 
finite  difference  zone  centers  adjacent  to  the  fault,  which 
are  0.05  km  from  the  fault  plane. 

First,  consider  the  stress  for  the  elastic  problem 
shown  in  Figure  15.  Initially,  the  stress  at  a given  point 
is  at  the  prestress  level  (1  kbar) . Prior  to  the  rupture 
front  arrival  at  a given  location,  the  stress  increases  above 
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the  prestress  level.  This  stress  concentration  ahead  of  the 
rupture  front  is  a general  characteristic  of  elas todynamic 
cracks  with  subsonic  rupture  propagation.  We  note  the  amplifi- 
cation of  this  stress  concentration  with  increasing  distance 
in  the  direction  of  rupture.  At  the  rupture  arrival  time, 
the  stress  drops  abruptly,  overshooting  and  then  settling  at 
the  prescribed  frictional  stress  of  820  bars.  The  overshoot 
results  from  the  fact  that  the  observation  points  are  slightly 
removed  from  the  fault  plane  itself.  The  same  phenomenon  is 
present  in  Richard's  (1976)  analysis  of  the  self-similar  ex- 
panding crack.  The  stress  then  remains  at  the  frictional 


stress  level  until  the  nearby  part  of  the  fault  plane  heals. 
Then  the  stress  relaxes  to  a value  less  than  that  of  kinetic 
friction.  The  fault  overshoots  the  static  equilibrium  value 
of  slip  (see  Section  3.5.3).  This  healing  wave  shows  up 
clearly  in  Figure  15;  it  propagates  toward  the  hypocenter  from 
the  periphery  of  the  fault.  The  last  phase  evident  in  the 
figure  propagates  outward  from  the  hypocenter,  and  this  phase 
corresponds  to  the  shear  wave  associated  with  the  final  ar- 
rest of  slip  at  the  center  of  the  fault. 

Now  consider  the  stresses  for  the  plastic  problem  in 
Figure  15.  Until  healing  occurs  at  a given  point,  the  stress 
history  is  unchanged  from  the  elastic  case  except  that  the 
stress  concentration  ahead  of  the  rupture  front  has  been  elim- 
inated. After  healing,  the  residual  stress  is  nearly  the 
same  as  in  the  elastic  case. 

Figure  16  displays  the  stress  near  the  fault  plane 

as  a fvmction  of  azimuth  4> . Again,  the  stress  is  zone 
centered,  so  it  is  actually  evaluated  at  points  0.05  km  dis- 
placed from  the  fault  plane.  The  five  points  displayed  are 
at  approximately  the  same  distance,  1.15  km,  from  the  center 
of  the  fault.  The  stress  concentration  preceding  rupture  is 
nearly  zero  (actually  slightly  negative)  at  <|)  = 0°,  and  in- 
creases smoothly  to  a maximum  of  (j)  = 90°.  This  pattern  can 
be  compared  to  Figure  8 of  Richards  (1976)  . The  stress 
histories  are  very  similar,  although  his  results  are  for  an 
ellipitcal  fault  in  which  rupture  -velocity  varies  with  azi- 
muth from  0.923  to  3/  whereas  our  numerical  solution  is  for 
circular  rupture  propagation  at  rupture  velocity  0.93. 

Figure  17  shows  the  stress  component  along  the 

fault  plane.  This  component  of  stress  is  not  relieved  by 
plastic  flow,  and  the  concentration  of  a ahead  of  the  rup- 
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ture  is  essentially  identical  in  both  the  plastic  and  the 
elastic  cases. 


I I 

1 second 

Figure  16.  Time  history  of  adjacent  to  the  fault  plane 

for  several  azimuthal  angles  41.  Each  observation 
point  is  at  a distance  of  0.05  km  from  the  fault 
plane,  eind  a distance  1.15  km  from  the  center  of 
the  fault. 


I 

I 


3.6.2  Particle  Velocity  on  the  Fault 

Figure  18  shows  the  particle  velocity  obtained  on  the 
fault  plane  at  increasing  distances  from  the  hypocenter  along 
a radial  line.  The  solid  curves  are  for  the  elastic  case, 
the  dashed  curves  are  for  the  plastic  case. 

It  is  evident  from  Figure  18  that  the  initial  velocity 

is  strongly  peaked  and  the  peak  value  increases  with  hypo- 

central  distance.  We  can  understand  this  characteristic  of 

the  velocity  curves  by  means  of  Kostrov's  analytic  solution. 

Equation  (23).  This  expression  implies  a velocity  propor- 
2 22  -1/2 

tional  to  t(t  - r /V_)  , which  is  singular  at  the  rupture 

K 

arrival  time  (except  at  the  hypocenter  r = 0) . Since  a 
finite  difference  grid  is  essentially  a low-pass  filter,  we 
expect  the  peak  velocities  in  Figure  18  to  be  a filtered  ver- 
sion of  the  analytic  solution.  To  approximate  the  filtering 
effect  of  the  grid  on  peak  velocity,  we  find  v,the  average  of 
the  velocity  over  a short  time  T following  the  rupture  ar- 
rival: 


C ^ ^ 
u T 


r2x-l/2 

h] 


= C — 6 t"^^^ 


which,  for  r/V  >>  T/2  is  proportional  to  r 
R 


That  is. 


we  expect  peaik  velocity  to  increase  as  the  square  root  of 
distance  in  the  direction  of  rupture.  The  peak  velocities 
in  Figure  19  are  in  good  agreement  with  Equation  (26)  if  T 
is  assumed  to  be  approximately  0.05  seconds  (ten  times  steps) 
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It  is  evident  from  Figure  18  that  the  initial  part  of 
the  particle  velocity  is  unaffected  by  the  introduction  of 
plasticity,  within  the  resolution  of  the  finite  difference 
calculation.  A large  velocity  peak  occurs  at  the  crack  tip, 
even  in  the  absence  of  the  strong  stress  concentration  asso- 
ciated with  the  purely  elastic  problem.  The  plastic  and 
elastic  solutions  are  indistinguishable  until  the  arrival  of 
the  stopping  phase. 

It  is  informative  to  compare  this  result  with  Ida's 
(1972,  1973)  analytic  study  of  propagating  cracks.  That 
study  was  an  elastodynamic  analysis,  but  finite  material 
strength  at  the  crack  tip  was  modeled  by  introducing  a dis- 
tribution of  "cohesive  stress"  across  the  crack  which  opposes 
the  slip.  This  approach  truncates  the  shear  stress  concen- 
tration at  the  crack  tip,  and  in  this  respect  simulates  the 
plastic  flow  model  used  here.  The  spatial  distribution  of  the 
cohesive  stress  must  be  assiimed  a priori,  however.  Ida  (1973) 
gives  an  approximate  expression  for  the  peak  particle  velocity: 


U 


peak 


a, 


(27) 


where  is  rupture  velocity,  is  the  excess  of  rock 
strength  over  sliding  friction.  The  D is  a constant  which 
was  found  to  be  about  one  from  numerical  experiments  with 
different  assumed  distributions  of  cohesive  stress.  For  our 
plastic  crack  model,  a is  180  bars,  V is  about  3 km/sec,  and 

c O R 

U is  about  3 x 10^  bars,  so  Equation  (27)  would  suggest  a 
value  of  cdaout  180  cm/sec  for  peak  velocity.  In  contrast,  the 
peak  velocity  for  the  finite  difference  solution  is  twice  this 
value  cind  is  a lower  limit  due  to  the  filtering  effect  of  the 
grid. 

The  discrepancy  is  due  to  the  fact  that  Ida  found  D 
to  be  of  order  1 by  considering  smoothly  varying  cohesive 
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. Slip  velocity  in  the  fault  plane.  Solid 
lines  are  the  elastic  case,  dashed  lines 
the  elastoplastic  case.  x,y  coordinates 
in  kilometers  are  shown  in  parenthesis. 


stress  distributions.  In  the  current  plastic  model, crack 
tip  stresses  are  limited  by  plastic  yielding  to  essentially 
a constant  (Figure  15) , until  the  rupture  arrival  time.  Then 
the  stress  abruptly  drops  to  the  kinetic  friction  level.  This 
discontinuous  behavior  of  the  stress  resembles  a special  case 
treated  by  Ida  (1972)  in  which  the  cohesive  stress  was  as- 
sumed to  behave  dis continuous ly . The  crack  tip  stress  was 
assumed  constant  over  a short  interval,  then  dropped  abruptly 
to  zero.  As  Equation  (22)  of  Ida  (1972)  indicates,  this 
special  case  yields  a singular  velocity  at  the  point  of  stress 
discontinuity,  in  spite  of  the  imposed  finiteness  of  stress. 

After  the  arrival  at  a given  point  of  edge  effects  due 
to  stopping  of  the  rupture  front,  the  fault  plane  velocities 
are  substantially  modified  by  the  plasticity.  As  Figure  18 
indicates,  yielding  at  the  crack  tip  smooths  the  stopping 
phase,  robbing  the  slip  function  of  high  frequencies  and  in- 
creasing the  long-period  content  of  the  slip  function.  The 
average  static  slip  on  the  fault  plane  is  increased  by  about 
11  percent  when  yield  is  permitted. 

3.7  FAR-FIELD  RESULTS 

In  this  section,  we  examine  the  radiated  seismic  field 
for  the  two  fault  models.  Far-field  displacement  spectra  and 
time  histories  are  given  for  P and  S waves.  The  far-field 
solutions  were  obtained  from  the  multipolar  expansions  (Sec- 
tion 4)  , using  terms  up  to  2.  = 8. 

3.7.1  Elastic  Case 

Figures  19  and  20  present  normalized  far-field  P and 
S wave  displacement  spectra  and  time  histories  for  the  elastic 
case.  Results  are  shown  at  10°  intervals  in  6;  Figure  19  is 
for  (p  = 90°,  and  Figure  20  is  for  <()  = 45°.  Solid  lines  are 
the  P wave  displacements,  dashed  lines  are  the  S wave  ones. 
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Normalizeci  far-field  P wave  (solid  curves)  and  S 
wave  (dashed  curves)  displacement  spectra  and  time 
histories.  Displacements  are  shown  at  10°  inter- 
vals in  9,  for  d)  = 90°. 
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Figure  20.  Normalized  far-field  wave  (solid  curves)  and  S 
wave  (dashed  curves)  displacement  spectra  and 
time  histories.  Displacements  are  shown  at  10 
intervals  in  0,  for  ()>  = 45°. 


The  curves  are  normalized  by  the  zero  frequency  spectral 
amplitude,  derived  from  the  average  slip  s.  For  P waves,  the 

For  S waves , we  can  in- 


normlization  is  yAs  (4TTpa^r) 


terpret  the  curves  either  as  the  9 component  of  displacement 

noinnalized  by  yAs  (4TTp6^r)  » or  as  the  <p  component  of 

3—1 

displacement  normalized  by  yAs  {4TTpg  r)  travel 

times  from  hypocenter  to  receiver  have  been  removed  from  the 
P and  S time  histories. 


Comparing  results  at  <j>  = 90°  with  those  at  ()>  = 45°, 
we  note  that  pulse  width  and  comer  frequency  have  practically 
no  dependence  on  <P ; at  higher  frequency  there  is  some  dif- 
ference in  spectral  and  time  domain  detail  between  the  two 
azimuths.  (Results  at  = 0°  are  virtually  identical  to 
those  at  (|>  = 90°,  and  are  now  shown.) 

Dependence  of  pulse  width  and  corner  frequency  on  0 
cind  on  wave  type  (P  or  S)  is  significant.  Our  observations 
concur  with  those  of  Madariaga  (1976) : 

(a)  S wave  corner  frequencies  are  smaller  than  P 
wave  corner  frequencies,  except  near  0=0°. 

(b)  Pulse  width  and  comer  frequency  are  governed 
by  the  travel  time  difference  between  stopping 
phases  from  the  near  and  far  edges  of  the  fault. 
Thus,  pulse  width  increases  with  0,  being 
greatest  for  observers  near  the  plane  of  the 
fault  and  smallest  for  observers  near  the  fault 
normal. 

(C)  P and  S wave  corner  frequencies,  for  0 > 30°, 

are  expressed  very  well  by  Madariaga's  Equation 
(24),  replacing  the  circular  fault  radius  with 
the  square  fault-width. 


3.7.2  Plastic  Case 

The  average  static  slip  for  the  plastic  fault  problem 
is  88  cm,  which  exceeds  the  average  slip  in  the  elastic  case 
by  11  percent.  This  results  from  the  smoothing  of  the  stopping 
caused  by  yielding  at  the  fault  edge.  Actually,  this  is 
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precisely  the  percentage  increase  that  would  be  predicted  by 
simple  scaling  of  the  elastic  problem  for  a fault  dimension 
occupying  the  entire  length  and  width  of  the  plastic  zone 
(the  plastic  zone  extended  0.2  km  beyond  edge  of  the  fault). 

The  seismic  moment  for  the  plastic  problem  cannot  be 
obtained  from  the  slip  function,  since  the  plastic  deforma- 
tion beyond  the  fault  edge  contributes  to  the  low-frequency 
spectrum  of  the  radiated  field.  Therefore,  we  derive  the 
seismic  moment  from  the  zero- frequency  limit  of  the  radiated 

field,  obtained  from  the  multipoles.  This  gives  a seismic 

2 4 

moment  of  3.15  x lo  dyne-cm,  which  is  38  percent  larger 
than  the  moment  for  the  elastic  case.  Scaling  of  the  elastic 
solution  as  suggested  in  the  last  paragraph  would  predict  a 
slightly  larger  increase  in  moment,  42  percent  instead  of  38 
percent. 

Figure  21  shows  the  effect  of  plastic  yield  on  the 
far-field  displacements.  Spectra  and  pulses  are  shown  at 
0 = 90®  for  three  values  of  9 : S wave  solutions  are  shown  at 
9=0®  and  9 = 90®,  and  P wave  solutions  are  shown  at  9 = 45® 
Dashed  curves  are  the  elastic  case,  solid  curves  the  plastic 
case.  In  each  case,  the  far-field  solution  is  the  sum  of 
multipolar  terms  up  to  il  = 8. 

The  main  influence  of  plastic  yielding  is  to  smooth 
the  stopping  phases,  with  the  result  that  tJle  low-frequency 

k 

part  of  the  spectrum  is  enhanced  at  the  expense  of  the  high 
frequency  part  of  the  spectrum.  Consider,  for  example,  the 
P wave  pulse  at  9 = 45®.  Two  stopping  phases,  corresponding 
to  rupture  arrival  at  the  near  and  far  edges  of  the  fault, 
respectively,  are  apparent  in  the  P wave  displacement  pulse 
for  the  elastic  problem.  These  appear  as  discontinuities  in 
slope  occurring  at  about  0.46  seconds  and  0.67  seconds.  The 
P wave  pulse  for  the  plastic  problem  coincides  with  that  for 
the  elastic  problem  until  the  arrival  of  the  first  stopping 
phase.  The  displacement  pulse  for  the  plastic  case  then 
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reverses  more  gradually,  overshooting  the  elastic  case;  the 
second  stopping  phase  is  almost  imperceptible  for  the  plastic 
case. 

Clearly,  the  increase  in  seismic  moment  is  the  conse- 
quence of  plastic  strain  induced  at  the  periphery  of  the  slip 
surface.  Relieving  stress  by  permitting  plastic  strain  is 
very  similar  to  relieving  stress  by  permitting  frictional 
sliding,  but  with  the  frictional  stress  approximately  equal 
to  Y//5  (which  in  this  case  equals  the  prestress  a^) . Out- 
side the  zone  in  which  yielding  was  permitted,  a static  shear 
stress  concentration  about  23  percent  in  excess  of  the  pre- 
stress developed.  Thus,  if  a larger  plastic  zone  had  been 
specified,  the  seismic  moment  would  have  been  even  greater, 
as  plastic  strains  extended  outward  to  eradicate  the  stress 
concentration.  On  the  other  hand,  had  a somewhat  higher 
yield  strength  been  specified,  the  results  of  the  plastic 
problem  would  have  approached  those  of  the  elastic  case.  If 
the  yield  parameter  Y//3  had  exceeded  the  tectonic  stress  by 
1.44  (a,j,  - a^)  (that  is,  if  Y had  been  26  percent  larger)  , 
no  yielding  would  have  occurred,  and  the  two  solutions  would 
have  been  indistinguishable. 

3.8  SUMMARY 

Two  three-dimensional  finite  difference  calculations, 
simulating  faulting  in  a uniform  whole  space,  have  been  per- 
formed. The  first  calculation  was  for  a linearly  elastic 
medium;  the  second  treated  an  elastic-plastic  medium.  Ac- 
curacy of  the  finite  difference  method  has  been  verified  by 
comparing  the  results  to  analytical  and  numerical  solutions 
i of  crack  problems.  The  far-field  displacements  associated 

I with  the  finite  difference  simulations  have  been  obtained 

I using  an  elastodynamic  representation  theorem.  This  integral 

, representation  has  been  applied  in  two  different  forms  and 

the  far-field  results  compared. 
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An  important  objective  was  to  evaluate  the  near-field 
and  far-field  effect  of  non-linear  deformation  in  the  fault 
zone.  It  was  found  that  the  initial  portion  of  the  slip 
time-function  is  unaffected  by  the  admission  of  a simple  form 
of  plasticity.  The  large  velocity  peak  at  the  crack  tip, 
characteristic  of  elastic  crack  problems,  persisted  in  the 
inelastic  case.  The  stopping  phase  was  modified  somewhat, 
however.  Yielding  at  the  edge  of  the  fault  resulted  in  less 
abrupt  stopping,  reducing  the  high-frequency  content  of  both 
the  slip  function  and  the  far-field  displacements,  and  in- 
creasing the  average  slip  by  11  percent.  Accumulation  of 
plastic  strain  beyond  the  fault  edge  resulted  in  a 38  percent 
higher  moment  thaui  in  the  elastic  case. 

These  effects  of  plasticity  depend  upon  our  choice  of 
the  magnitude  of  prestress,  the  yield  strength,  the  fric- 
tional stress,  and  the  width  of  the  non-linear  zone.  The 
problem  treated  was  an  extreme  case  in  the  sense  that  the 
prestress  level  everywhere  equaled  the  strength  of  the  medium 
(Y//3) . With  a moderate  increase  in  the  yield  strength  (26 
percent),  or  a similar  decrease  in  the  prestress,  there  would 
have  been  no  yielding,  and  the  solution  would  have  been 
identical  to  that  of  the  elastic  problem.  On  the  other  hand, 
we  arbitrarily  limited  the  extent  of  the  plastic  zone.  Out- 
side tne  plastic  zone,  a static  shear  stress  concentration 
persisted  which  was  nearly  as  large  as  that  of  the  elastic 
case;  had  the  plastic  zone  been  larger,  more  plastic  strain  j| 

would  have  occurred,  resulting  in  an  even  larger  seismic  'j 

moment.  ' i 
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IV.  COMPARISON  OF  ANALYTICAL  AND  FINITE  DIFFERENCE 

MODELS 

4.1  INTRODUCTION 

In  Section  II  we  developed  a model  for  the  Pocatello 
Valley  earthquake  which  is  expressed  in  terms  of  the  analyti- 
cal relaxation  model  of  Archambeau  (1968)  and  Minster  (1973) . 
This  analytical  model  is  convenient  to  use  because  it  is  ex- 
pressed in  terms  of  relatively  few  parameters.  However, 
interpretation  of  these  parameters  in  terms  of  earthquake 
physics  is  quite  controversial  because  of  the  geometry  in 
which  the  source  model  is  formulated.  The  Archambeau/Minster 
source  is  a spherical  volume  of  reduced  shear  strength  which 
grows  in  a pure  shear  prestress  field.  The  growth  is  asym- 
metric with  a point  on  the  boundary  being  fixed  as  shown  below. 
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This  model  does  include  directional! zed  rupture,  finite  source 
dimension  and  finite  rupture  velocity  effects,  but  in  a geometry 
unlike  conventional  ideas  of  earthquake  geometries. 

The  ILLIAC  calculations  of  earthquaike  faulting  are  close 
to  the  most  detailed  and  realistic  models  of  earthquake  fault- 
ing that  are  currently  available.  The  disadvantage  of  using 
such  calculations  to  model  a particular  earthquake  is , of 
course,  the  expense  of  varying  the  parameters.  The  parameters 
for  our  two  calculations  are  not  entirely  fixed,  however,  be- 
cause the  fault  dimension  and  stress  drop  can  be  scaled. 
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Taken  together,  the  analytic  models  and  the  ILLIAC 
calculations  can  be  used  to  define  the  earthquake  source  with 
considerable  confidence.  Questions  about  the  physical  signif- 
icance of  parameters  of  the  analytical  model  Ccui  be  resolved 
by  comparing  to  the  numerical  model  results. 

In  this  section  we  compare  the  analytical  model  for  the 
Pocatello  Valley  earthquake  from  Section  II  with  the  finite 
difference  model.  The  implications  of  this  comparison  for 
the  Pocatello  Valley  event  and  for  western  U.S.  earthquakes 
in  general  will  also  be  discussed. 


4.2  SCALING  THE  FINITE  DIFFERENCE  MODEL 

In  Section  3.3  we  listed  the  parameters  for  the  two 
ILLIAC  earthquake  calculations.  For  the  elastic  case  all 
parameters  could  be  nondimens ionali zed  and  the  results  could 
have  been  presented  in  that  form.  The  elas toplastic  fault 
zone  calculation  is  not  so  easily  scaled  and  so  we  presented 
all  our  results  in  terms  of  specific  parameters. 

For  the  elastic  problem,  the  fundamental  parameters  are 
the  aspect  ratio  of  the  fault  (in  this  case,  it  is  square) , 
Poisson’s  ratio  (0.25),  and  the  ratio  Vj^/6  (in  this  case  0.9). 
If  these  dimensionless  parameters  are  held  fixed,  the  ordinate 
and  abiL-cissa  of  the  far-field  displacement  spectrum  can  be 
scaled  with  fault  dimension  a,  stress  drop  Aa , shear  modulus 
p,  shear  speed  8,  and  hypocentral  distance  r.  A nondimensional 
frequency  f and  nondimensional  moment  have  the  form 


f ' 


fa 

B 


M' 


o 


M 

o 

Aoa^ 


(28) 


Thus,  the  far-field  spectral  displacement  scales  as  a^Ao/ySr 
and  frequency  scales  as  3 /a. 
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The  wave  velocity  parameters  a and  6 and  the  shear 
modulus  p are  appropriate  for  the  source  region  of  the 
Pocatello  Valley  earthquake  and,  in  fact,  for  a majority  of 
the  shallow  western  U.S.  earthquakes.  Therefore,  we  leave 
these  parameters  fixed.  The  rupture  velocity  is  then  fixed 
to  be  V = 3.08  km/sec  or  0.9B.  The  parameters  we  are  free 
to  scale  are  the  fault  dimension,  a,  and  the  stress  drop.  La. 

In  order  to  compare  to  data,  we  might  wish,  for  ex- 
cimple,  to  scale  the  comer  frequency  without  modifying  the 
long  period  spectral  level  (the  seismic  moment).  To  change 
the  comer  frequency  by  a factor  q,  we  would  replace  the 
source  dimension  a by  a/q  and  replace  the  stress  drop  La  by 

3 

Aa/q  . Similarly,  to  multiply  the  moment  by  a factor  of  q 
without  modifying  the  corner  frequency,  we  would  replace  La 
with  qAa , and  leave  the  source  dimension  a unchanged. 

The  elastoplas tic  source  is  not  so  easily  scaled. 
Stress  drop  scaling  can  only  be  approximate.  Zones  that  al- 
most yield  for  one  stress  drop  will  yield  for  a higher  stress 
drop.  On  the  other  hand,  lowering  the  stress  drop  will  re- 
move yielding  from  some  zones.  It  is  not  known  how  large  an 
effect  this  is.  However,  we  note  that  the  stress  drop  scal- 
ing is  generally  by  a large  factor  since  it  is  proportional 
to  the  cube  of  the  length  scaling. 

What  about  the  length  scaling  of  the  elastoplas tic 
case?  We  Ccin  scale  exactly  with  length  according  to  (28)  as 
long  as  the  cuiswers  are  independent  of  the  size  of  the  plas- 
tic zone.  That  is,  if  the  stresses  outside  the  plastic  zone 
are  everywhere  less  than  the  yield  condition.  However,  as 
discussed  in  Section  3.7.2,  the  displacement  field  for  this 
elas toplastic  calculation  is  dependent  on  the  dimensions  of 
the  plastic  zone.  The  quantitative  extent  of  this  dependence 
is,  however,  not  known. 


There  is  some  difficulty  in  obtaining  a suitable  defi- 
nition of  Aa.  The  quantity  - a the  difference  between 
tectonic  and  frictional  shear  stress,  is  a well-defined 
physical  par2uneter  in  our  fault  model.  However,  we  would 
prefer  to  define  Aa  to  represent  a static  stress  drop,  since 
it  would  then  be  comparable  to  the  parameter  appearing  in  the 
analytic  fault  model.  The  static  stress  drop  is  a well-de- 
fined quantity  at  a given  location  on  the  fault  plane,  namely 
the  difference  between  the  initial  and  final  shear  stress. 
However,  this  quantity  is  not  a constant,  but  varies  with 
position  on  the  fault  (see  Figure  16) . Nor  is  the  spatial 
average  of  the  static  stress  drop  necessarily  a meaningful 
quantity,  since  a large  static  stress  concentration  over  a 
small  area  can  greatly  affect  this  average,  with  minimal  ef- 
fect on  the  radiated  field.  Instead,  we  shall  define  Aa  to 
be  1.14  (a^  - a^) ; this  choice  of  Aa  gives  the  correct  value 
of  seismic  moment  when  the  average  slip  is  estimated  from 
Equation  (24)  and  then  the  moment  obtained  from  Equation  (25) 
(see  Section  3.5.3).  With  this  definition,  Aa  for  the  elas- 
tic finite  difference  model  is  205  bars.  We  note  that  the 
actual  static  stress  drop  on  the  fault  varies  fairly 
smoothly;  it  is  1.29  (a^  - a^)  at  the  center,  decreasing  to 
approximately  1.1  (a^  - a^)  at  a distance  of  2a/3  from  the 
center,  with  little  azimuthal  variation. 

4. 3 ELASTIC  FINITE  DIFFERENCE  SOURCE  COMPARED  TO  OBSER- 
VATIONS OF  THE  POCATELLO  VALLEY  EARTHQUAKE 

Bache  and  Harkrider  (1976)  describe  a technique  for 
representing  the  output  of  finite  difference  source  calcula- 
tions like  those  in  Section  III  in  terms  of  a multipolar  ex- 
pansion. This  technique  is  summarized  in  Section  3.4.  Once 
we  have  computed  the  expansion  coefficients,  the  multipole 
coefficients,  in  the  source  coordinate  system,  the  coordinate 
rotation  formulae  given  by  Minster  (1976)  can  be  used  to  ro- 
tate the  source  to  any  desired  orientation. 
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The  Archcunbeau/Minster  source  used  in  Section  II  to 
develop  the  analytical  model  for  the  Pocatello  Valley  earth- 
quake is  expressed  directly  in  terms  of  multipole  coefficients 
Thus , both  sources  can  be  represented  in  the  same  way  and 
direct  comparison  is  facilitated.  Synthetic  seismograms  are 
computed  using  the  methods  and  earth  models  described  in 
Section  2.4. 

The  long  period  body  wave  seismogreims  for  the  finite 
difference  sources  are  essentially  identical  to  those  shown 
in  Figure  7.  For  these  data  the  only  important  source  pareim- 
eters  are  the  moment,  focal  depth  and  orientation. 

Recall  from  Section  3.5.3  that  the  moment  for  the 

elastic  finite  difference  calculation  with  the  parameters 

24 

given  in  Section  3.3  is  2.28  x 10  dyne-cm.  The  moment  for 

the  Pocatello  Valley  earthquake  inferred  from  the  long  period 

2 4 

body  waves  is  7.0  x lo  dyne-cm.  Thus,  we  need  to  scale  the 
calculation  to  a moment  that  is  larger  by  a factor  of  3.1. 

The  most  straightforward  way  to  scale  the  elastic 
finite  difference  source  to  the  moment  of  the  Pocatello  Valley 

I 

event  is  to  increase  the  stress  drop  from  205  bars  to  630 
bars.  Let  us  do  so  and  explore  the  consequences  for  the 
short  period  seismogreims. 

In  Figure  22  we  add  the  short  period  seismograms  for 
the  elastic  finite  difference  source  with  Aa  = 630  bars  to 
Figure  8 from  Section  II.  The  waveforms  and  frequency  con- 
tent are  rather  close  to  those  for  Model  I,  the  analytical 
source  with  a fault  length  of  3 km.  The  main  difference  is 
that  for  the  finite  difference  source  the  P phase  is  larger 
relative  to  the  pP  phase.  This  is  a consequence  of  the  bi- 
lateral nature  of  the  finite  difference  source,  as  we  shall 
see. 

The  Archambeau/Minster  analytical  model  can  be  modi- 
fied to  approximate  bilateral  faulting.  We  do  so  by 
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Figure  22.  Short  period  seismograms  from  the  elastic  finite  difference  source  are 
compared  to  the  observations  of  the  Pocatello  Valley  earthquake  and  the 
synthetic  records  from  two  analytical  source  models  (Figure  8) . 


superimposing  the  radiation  fields  from  two  equal  unilateral 
faults  rupturing  in  opposite  directions  in  the  same  prestrain 
field.  The  approximation  is  that  the  radiation  field  from 
each  unilateral  segment  does  not  include  effects  due  to  the 
presence  of  the  other  segment.  The  mathematical  consequences 
of  approximating  bilateral  faulting  in  this  way  are  quite 
elementary  and  are  worth  mentioning.  The  multipolar  expan- 
sion of  the  radiation  field,  Equation  (17)  in  Section  3.4.1, 
includes  terms  even  and  odd  in  H . The  dominant  term  at  long 
period  is  the  double-couple , 1=  2.  The  other  even  order 
terms  account  for  finite  source  effects  that  are  independent 
of  rupture  direction.  The  odd  order  terms  account  for  ef- 
fects due  to  unilateral  rupture  direction.  Thus,  a unilateral 
model  like  Model  I Ccin  be  made  bilateral  by  merely  dropping 
the  odd  order  terms  from  the  calculation  cind  doubling  the 
even  order  terms. 

Let  us  now  directly  compare  the  seismograms  from  the 

elastic  finite  difference  source  to  those  from  the  unilateral 

and  bilateral  versions  of  Model  I.  All  three  sources  are 
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scaled  to  the  same  moment  which  is  7.0  =<  10  dyne-cm.  The 
finite  difference  source  has  faulting  on  a 3 km  x 3 km  rec- 
tangular fault  plane.  The  fault  dimension  for  Model  I is 
L = 3 km,  where  L is  the  diameter  of  the  spherical  volume  in- 
to which  stress  is  released.  The  bilateral  version  of  Model 
I has  two  such  spherical  volumes.  The  planar  area  of  fault- 
ing (ignoring  the  out-of-plane  dimension)  and  stress  drop  for 
each  model  are  summarized  below: 

Planar  Area  Stress  Drop, 

of  Faulting  (km^)  Aa  (bars) 

Elastic  Finite 

Difference  Source  9 630 

Unilateral  Model  I 7.1  253 


Bilateral  Model  I 


14.1 


127 


A bilateral  Archambeau/Minster  model  with  a planar 

fault  area  of  9 km  would  include  two  sections,  each  with 

L = 2.4  km.  The  stress  drop  which  gives  the  same  Aa  as 

the  bilateral  Model  I is  250  bars.  A unilateral  Archambeau/ 

2 

Minster  model  with  a planar  area  of  9 km  has  L = 3.4  km  and 
Aa  = 176  bars.  These  stress  drops  are  factors  of  2.5  and  3.6 
times  smaller  than  the  stress  drop  associated  with  the  rec- 
tangular fault.  These  factors  compare  to  the  value  of  3.6 
which  scales  the  stress  drop  of  the  Archambeau/Minster  model 
to  that  of  a circular  fault  of  the  same  planar  area  and  same 
moment  (Minster  and  Suteau,  1976) . 

In  Figure  23  we  compare  the  five  short  period  tele- 
seisms from  the  elastic  finite  difference  calculation  to  those 
from  the  unilateral  and  bilateral  versions  of  Model  I.  The 
amplitude  (corrected  for  period  dependent  instrument  re- 
sponse) and  period  data  from  these  synthetic  seismograms  are 
summarized  in  Table  7.  The  bilateral  Model  I is  somewhat 
longer  period  than  the  elastic  finite  difference  source, 
which  is  expected  since  it  has  larger  source  dimension.  The 
closer  agreement  of  the  bilateral  Model  I with  the  finite 
difference  source  is  most  apparent  in  the  ratio  of  pP  to  P 
phases , though  this  is  more  clear  from  the  records  in  Figure 
23  thsn  from  the  numbers  in  the  table. 

Now  let  us  return  to  the  comparison  of  the  seismograms 
from  the  finite  difference  source  with  the  Pocatello  Valley 
earthquake  observations.  Figure  22.  In  Table  8 we  compare 
the  amplitudes  and  periods  in  the  same  format  used  for  Tcibles 
5 cuid  6 in  Section  II.  We  see  that  we  have  precisely  the  scime 
problem  we  had  with  Model  I in  Section  II.  That  is,  when  the 
finite  difference  source  is  scaled  to  have  the  correct  moment, 
the  short  period  seismograms  are  too  large  by  a factor  of 
four  or  five. 
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TABLE  7 


AMPLITUDES  AND  PERIODS  OF  THE  SEISMOGRAMS  OF  FIGURE  23 


P PHASE 

MAXIMUM 

PHASE 

STATIONS 

Amplitude 

(microns) 

Period 

Amplitude 

(microns) 

Period 

A„/A 

P max 

RES 

2.11 

1.1 

2.59 

1.0 

1.23 

EFD  Source 

Uni-Model  I 

1.80 

1.2 

2.33 

1.0 

1.29 

Bi-Model  I 

1.84 

1.3 

1.96 

1.2 

1.07 

DAG 

EFD  Source 

1.38 

1.1 

2.38 

1.2 

1.73 

Uni -Model  I 

1.14 

1.2 

2.27 

1.2 

1.99 

Bi -Model  I 

1.19 

1.2 

o T 

M • X ^ 

1.3 

1.78 

KTG 

EFD  Source 

1.29 

1.1 

2.42 

1.2 

1.88 

Uni -Model  I 

1.07 

1.2 

2.37 

1.2 

2.22 

Bi -Model  I 

1.11 

1.2 

2.24 

1.3 

BOG 

EFD  Source 

0.76 

1.1 

1.75 

1.3 

2.30 

Uni -Model  I 

0.72 

1.2 

1.72 

1.2 

2.39 

Bi -Model  I 

0.69 

1.2 

1.53 

1.4 

2.21 

PTO 

EFD  Source 

0.67 

1.1 

1.96 

1.3 

2.93 

Uni-Model  I 

0.72 

1.3 

2.02 

1.3 

2.81 

Bi -Model  I 

0.61 

1.2 

1.88 

1.4 

3.08 

*Elastic  Finite  Difference  Source  Model 


TABLE  8 

AMPLITUDE  AND  PERIOD  COMPARISONS  FOR  THE  ELASTIC  FINITE 
DIFFERENCE  SOURCE  SEISMOGRAMS  IN  FIGURE  22 


P PHASE 

MAXIMUM  AMPLITUDE 

Station 

Amplitude  Ratio 
Aq/As 

Periods 

Tq/Ts 

Amplitude  Ratio 
Aq/As 

Periods 

To/Ts 

RES 

0.06 

1.2/1. 1 

0.15 

1. 3/1.0 

DAG 

0.07 

1. 1/1.1 

0.20 

1. 3/1.2 

KTG 

0.15 

1. 1/1.1 

0.33 

1.2/1. 3 

BOG 

0.35 

1. 4/1.1 

0.33 

1.2/1. 3 

PTO 

0.20 

1. 4/1.1 

0.18 

1. 2/1.0 

Logarithmic 

Mean 

0.13 

0.22 

Standard 

Deviation 

106% 

41% 

Ccin  the  finite  difference  source  be  scaled  to  have 
less  high  frequency  energy  while  retaining  the  same  moment? 
The  appropriate  scaling  relations  are  described  in  Section 
4.2.  Let  us  consider  three  earthqucike  sources: 


1. 

a = 

1.5  km, 

Aa  = 6 30 

bars ; 

2. 

a = 

2.25  km. 

Aa  = 187 

bars ; 

3. 

a = 

3 km. 

Aa  = 79 

bars . 

The  first  is  the  elastic  finite  difference  source  with  the 
original  fault  dimension  (a  is  defined  in  Figure  12) . The 
far-field  displacement  spectra  for  this  source  appear  in 
Figure  20.  The  second  and  third  sources  have  the  same 
spectra  except  that  the  frequency  axis  is  scaled  by  1.5  and 
2.0,  respectively.  That  is,  the  corner  frequency  moves  to 
lower  frequency.  All  three  sources  have  the  same  moment. 

In  Figure  24  we  compare  the  seismograms  at  Station  ■ 

KTG  for  these  three  source  models.  The  amplitude  and  period  ' 

data  from  these  seismograms  are  compared  in  Table  9.  We  see  ' 

that  the  effect  of  moving  the  corner  frequency  to  lower  fre-  I 

quency  is  quite  large  on  the  period.  However,  the  short  ; 

period  amplitudes  do  not  get  much  smaller!  i 

I 

In  Figure  25  we  compare  the  seismograms  at  all  five  | 

static  .is  for  the  first  and  third  models.  The  aimplitude  and 

period  data  are  compared  in  Table  10.  At  all  stations  we  j 

see  that  the  measured  periods  are  0.4-0. 8 seconds  longer  for 

the  larger  source  dimension  fault.  However,  the  amplitudes 

(after  correction  for  frequency  dependent  instrument  re-  ' 

sponse)  are  not  much  different  for  the  two  sources.  | 

( 

The  results  of  this  scaling  of  the  elastic  finite  dif-  ! 

ference  source  support  our  conclusions  about  earthquake  fault- 
ing reached  earlier  for  the  San  Fernando  (Bache  and  Barker, 

1978)  and  Pocatello  Valley  (Section  II)  earthquakes  follow- 
ing our  analysis  with  the  analytical  Archeimbeau/Minster  model. 
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! Figure  24.  Synthetic  seismograms  are  compared  at  Station  KTG 

for  the  elastic  finite  difference  source  scaled  to 
three  source  dimensions  with  the  source  moment 
held  constant. 


TABLE  9 


AMPLITUDE  AND  PERIOD  DATA  FOR  SEISMOGRAMS  AT  KTG  FOR  THREE 
SCALED  VERSIONS  OF  THE  FINITE  DIFFERENCE  SOURCE  (FIGURE  24) 


P PHASE 

MAXIMUM 

PHASE 

Source 

Amplitude 

(microns) 

Period 

(seconds^ 

Amplitude 

(microns) 

Period 

(seconds) 

a = 

1.5,  Aa=  630 

1.29 

1.15 

2.42 

1.21 

a = 

2.25, Aa=  187 

1.04 

1.30 

2.21 

1.41 

a = 

3,  Aa=  79 

0.95 

1.56 

2.55 

1.79 
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AMPLITUDE  AND  PERIOD  DATA  FOR  TWO  SCALED  VERSIONS  OF  THE 
ELASTIC  FINITE  DIFFERENCE  SOURCE  (FIGURE  25) 


That  is,  we  cannot  simultaneously  match  the  long  and  short 
period  data  with  a single  rupture  velocity /single  stress  drop 
source  model.  In  the  following  section  we  will  see  what  ef- 
fect the  inclusion  of  elastoplastic  material  behavior  will 
have. 


4.4  ELASTOPLASTIC  FINITE  DIFFERENCE  SOURCE  COMPARED  TO 

OBSERVATIONS  OF  THE  POCATELLO  VALLEY  EARTHQUAKE 

In  Section  3.7.2  we  compared  the  far-field  displacement 
spectra  from  the  elastic  and  elastoplastic  finite  difference 
sources.  The  differences  were  that  the  elastoplastic  source 
has  larger  moment  and  lower  comer  frequency.  These  dif- 
ferences make  the  elastoplastic  source  look  much  like  the 
elastic  source  scaled  to  a larger  source  dimension,  at  least 
as  far  as  the  amplitude  spectra  are  concerned.  Let  us  now 
study  the  seismograms  from  the  elastoplastic  source. 

The  moment  of  the  elastoplastic  source  as  computed  is 
24 

3.15  X 10  dyne-cm.  This  is  smaller  than  the  moment  of  the 
Pocatello  Valley  earthquake  by  a factor  of  2.2.  We  pointed 
out  in  Section  4.2  that  the  elastoplastic  calculations  cam 
only  be  scaled  approximately.  But  since  comparison  is  easier 
if  all  the  synthetic  seismograms  have  the  saune  moment,  let 
us  scale  the  stress  drop  to  451  bars,  i.e. , by  a factor  of 
2.2.  The  scaled  source  spectra  are  probably  not  much  dif- 
ferent than  we  would  have  gotten  if  we  had  repeated  the  cal- 
culation with  this  new  stress  drop. 

In  Figure  26  we  compare  the  observations  with  the 
synthetic  seismograuns  from  the  elastic  and  elastoplastic 
source  calculations  with  the  original  source  dimension  (a  = 

1.5  km).  For  the  elastoplastic  source  the  amplitudes  and 
periods  are  compared  to  the  observations  in  Taible  10  in  the 
usual  format. 

From  the  seismograms  in  Figure  26  and  the  data  in 
Table  11,  we  see  that  the  addition  of  elastoplastic! ty  to  the 
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Figure  26.  Seismograms  from  the  two  finite  difference  source  calculations  (scaled  by 
stress  drop  to  the  correct  moment)  are  compared  to  observations  of  the 
Pocatello  Valley  earthquake. 


AMPLITUDE  AND  PERIOD  COMPARISONS  FOR  THE  ELASTOPLASTIC  FINITE 
DIFFERENCE  SOURCE  SEISMOGRAMS  IN  FIGURE  26 


P 

PHASE 

MAXIMUM  AMPLITUDE 

Station 

Amplitude 

Ao/Ag 

Ratio 

Periods 

To/Tg 

Amplitude  Ratio 

Aq/ Ag 

Periods 

To/Tg 

RES 

0.07 

1.2/1. 2 

0.16 

1. 3/1.1 

DAG 

0.07 

1.1/1. 2 

0.21 

1.3/1. 3 

KTG 

0.16 

1.2/1. 2 

0 .31 

1.3/1. 3 

BOG 

0.36 

1. 4/1.1 

0 .36 

1.2/1. 3 

PTO 

0 .21 

1. 4/1.1 

0 .20 

1. 2/1.0 

Logarithmic 

Mean  0 . 14 

0 .24 

Standard 

Daviation  104% 

39% 

model  causes  the  dominant  periods  to  be  a bit  larger,  about 
0.1  seconds  in  most  cases.  With  the  moment  at  the  right 
level  the  synthetic  short  period  seismograms  are  still  too 
large  by  a factor  of  four  or  more. 

4.5  CONCLUSIONS 

The  comparison  of  the  analytical  and  finite  difference 
models  confirms  our  most  important  conclusion  in  Section  II. 
The  long  and  short  perior  data  cannot  be  simultaneously  fit 
with  a single  rupture  velocity,  single  stress  drop  model. 
Variable  rupture  velocity,  variable  stress  drop  effects  are 
important. 

In  this  section  we  have  studied  two  finite  difference 
earthquake  simulations.  Since  these  are  for  simple  events, 
they  are  not  able  to  match  the  long  and  short  period  data 
simultaneously.  When  we  scale  them  to  the  right  moment,  the 
short  period  radiation  is  too  large  by  a factor  of  four  or 
more. 

Comparing  seismograms  from  the  analytical  and  finite 
difference  models  we  see  the  importance  of  source  directivity 
for  matching  the  data.  The  finite  difference  sources  are 
bilateral  and  give  poorer  agreement  with  the  data  than  the 
unilateral  analytical  source. 

Another  important  and  interesting  facet  of  the  com- 
parison is  that  it  gives  an  opportunity  to  study  the  physical 
meaning  of  the  parameters  of  the  Archambeau/Minster  model 
which  is  cast  in  terms  of  an  unrealistic  spherical  geometry. 
The  radiation  pattern  of  P and  S waves  from  the  analytical 
model  is  not  the  same  as  from  a planar  carck,  but  the  dif- 
ferences do  not  appear  to  be  too  important. 

The  physical  meaning  of  the  important  parameter  stress 
drop  is  made  more  clear  after  comparing  the  cinalytical  and 


finite  difference  models.  We  can  scale  the  analytical  model 
to  have  the  scune  planar  area  as  the  finite  difference  fault 
plane.  For  the  two  models  to  have  the  same  area,  the  stress 
drop  for  the  analytical  source  is  a factor  of  2.5  (bilateral) 
or  3.6  (unilateral)  times  smaller  than  for  the  finite  dif- 
ference source.  Thus,  stress  drops  are  underestimated  by 
eibout  this  cimount  if  they  are  based  on  the  Archcimbeau/ 

Minster  model.  This  confirms  the  stress  drop  estimate  for  the 
initial  faulting  of  the  Pocatello  Valley  event  from  Model  II 
given  at  the  end  of  Section  II;  that  is,  162  bars. 
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